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linear finite  difference  simulations  of  stick-slip  earthquake  faulting 
and  an  explosion  in  an  axesymmetric  finite  length  tunnel. 

In  Section  IV,  both  generalized  ray  theory  and  matrix  techniques 
were  used  to  calculate  synthetic  body  wave  seismograms  for  a point  source 
in  a layered  elastic  medium.  Both  methods  are  found  to  be  equivalent 
and  complementary  for  various  model  calculations.  Ray  theory  has  the 
advantage  of  giving  a direct  physical  interpretation  to  observed  phases. 
Matrix  techniques,  while  exact  in  terms  of  including  all  multipoles, 
give  little  insight  into  the  model,  but  serves  as  a check  on  the 
dominant  phases  needed  in  the  ray  calculations. 


In  Section  V,  a data  set  of  forty-one  moderate  to  large  earthquakes 
were  used  to  derive  scaling  rules  for  kinimatic  fault  models.  ^If 
effective  stress  and  static  stress  drop  are  equal,  then  fault  rise  time 
t,  and  fault  area,  S,  are  related  by  t = 16/T!>/(7tt  ^8) . Fault  length 

and  width  are  empirically  related  by  L = 2W.  These  scaling  laws  combine 
to  give  width  and  rise-time  in  terms  of  fault  length.  Finally  averaging 
reported  rupture  velocities  yields  V„  = .72B. 
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I. 


SUMMARY 


Previous  to  14  March  1975,  the  principal  investigator  for  this  contract 
was  Don  L.  Anderson,  with  co-investigators  C.  B.  Archambeau,  D.  G.  Harkrider, 
D.  V.  Helmberger  and  H.  Kanamori.  The  title  during  this  period  was  "Seismic 
Phenomena  Connected  with  Earthquakes  and  Explosions." 

The  major  accomplishments  during  this  period  can  be  broadly  categorized 
under  the  headings  of  (1)  Applications  of  1. inear  Inversion  Theory, (2)  Seismic 
Source  Theory,  and  (3)  Synthetic  Body  and  Surface  Wave  Modeling  of  Observed 
Seismic  Events. 


(1)  Application  of  Linear  Inversion  Theory 

Linear  inversion  theory  was  applied  to  body  wave  travel  times  and  free 
oscillation  data  to  obtain  the  velocity  and  density  structure  of  a monopole 
earth.  The  theory  was  adapted  to  find  the  set  of  source  parameters  which 
optimally  fit  the  variance  weighted  data  set  of  spatial  and  frequency 
distributed  path-corrected  Rayleigh  waves  from  the  San  Fernando  earthquake. 

The  interesting  side  result  from  this  technique  is  the  a priori  determination 
of  what  azimuth  distribution  of  observations  is  most  beneficial  in  determining 
the  source  parameters.  Of  the  18  WWSN  stations  used  for  this  event,  two 
stations  contain  a total  of  30%  of  the  total  information  in  the  data  set. 

(2)  Seismic  Source  Theory 

Two  different  approaches  for  investigating  the  characteristics  of  the 
seismic  source,  in  particular  the  determination  of  corner  frequency  for  P 
and  S waves,  were  developed.  The  first  uses  theoretical  volume  source  models 
of  tectonic  release  of  the  Archambeau  tvpe  and  obtains  numerically  their 
P-  and  S-wave  spectra.  The  spectra  for  these  particular  models  show  a definite 
difference  in  the  P-waves  and  S-wave  corner  frequencies.  This  approach  yields 
theoretical  Mg  vs  mp  curves  which  can  be  discussed  in  terms  of  scaling  by  the 
various  source  parameters.  The  second  approach  is  to  obtain  the  earthquake 
P-wave  source  history  from  comparison  of  theoretical  and  observed  seismograms. 
This  history  is  then  used  as  the  S-wave  source  history  to  obtain  theoretical 
seismograms  which  are  compared  with  S-wave  observations  of  the  same  events. 

The  question  of  the  validity  of  using  far-field  theory  in  order  to 
determine  source  parameters  such  as  moment  and  corner  frequencies  at  ranges 
where  this  assumption  is  questionable  was  evaluated.  The  displacement  fields 
were  numerically  investigated  at  a variety  of  ranges  and  source  histories  in 
order  to  establish  criteria  for  estimating  the  minimum  far  field  range  in 
both  the  time  and  frequency  domain. 

The  effect  of  finite  source  regions  on  various  scaling  laws  for  surface 
waves  from  underground  explosions  were  investigated  using  simple  theoretical 
source  models.  The  scaling  laws  for  simple  models  were  found  to  be  critically 
influenced  by  the  effective  radius  of  t lie  source  region.  The  strong  dependence 
on  radius  appears  to  he  related  to  the  perfect  symmetry  of  the  source  regions 
and  will  be  investigated  later. 
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3)  Synthetic  Body  and  Surface  Wave  Modeling  of  Observed  Seismic  Events 

Amplitudes  and  arrival  times  of  body  and  Rayleigh  waves  from  NTS 
events  were  used  to  determine  an  average  crustal  model  from  NTS  to  mid-Arizona. 
Transversely  polarized  shear  waves  (SH)  from  well  located  earthquakes  were  also 
used  to  determine  an  upper  mantle  shear  velocity  structure.  It  was  further 
demonstrated  that  observations  of  travel  times  for  vertically  polarized  shear 
waves  (SV),  the  most  common  technique,  can  lead  to  erroneous  shear  structure 
models  due  to  contamination  from  P-SV  interactions. 

In  order  to  study  body  waves  from  realistic  models  of  earthquakes, 
generalized  ray  theory  was  extended  to  a shear  dislocation  source. 

Displacements  at  the  surface  of  layered  halfspace  produced  by  a point 
source  dislocation  were  investigated.  Expressing  the  source  in  terms  of 
P,  SV  and  SH  displacement  potentials  allows  the  solution  to  be  expanded  in 
generalized  rays.  The  transient  response  for  each  ray  is  obtained  by  the 
Cagniard-deHoop  method.  First-motion  and  high-frequency  asymptotic 
approximations  of  the  exact  solutions  were  discussed. 

Using  a generalized  inverse  technique  WWSSN  long-period  P and  SH 
waveforms  were  analyzed  from  the  Koyna  earthquake.  The  effects  of  local 
plane-layered  earth  structure  near  an  imbedded  point  dislocation  source 
are  put  in  using  a modified  plane-wave  ray  theory  which  includes  the 
standard  reflection  and  transmission  coefficients  plus  source  corrections 
for  radiation  pattern  and  geometrical  spreading.  The  generalized  inverse 
compares  synthetic  seismograms  to  the  observed  in  the  time  domain  through 
the  use  of  a correlation  function.  Using  published  crustal  models  of  the 
Koyna  region  and  primarily  by  modelling  the  crustal  phases  P,  pP,  and  sP, 
the  first  25  seconds  of  the  long-period  waveforms  are  synthesized  for 
17  stations  and  a focal  mechanism  obtained  for  the  Koyna  earthquake  which 
is  significantly  different  from  previous  mechanisms. 

The  technique  has  also  been  successfully  applied  to  the  Borrego 
Mountain  earthquake  of  April  9,  1968.  Synthetic  seismograms  computed  from 
the  resulting  model  match  in  close  detail  the  first  25  seconds  of  long 
period  seismograms  from  a wide  range  of  azimuths.  The  main  shock  source 
time  function  was  determined  by  a new  simultaneous  short  period-long  period 
deconvolution  technique  as  well  as  by  the  inversion  technique.  The  duration 
and  shape  of  this  time  function  indicate  that  most  of  the  body  wave  energy 
was  radiated  from  a surface  with  effective  radius  of  only  8 km.  This  is 
much  smaller  than  the  total  surface  rupture  length  or  the  length  of  the^,. 
aftershock  zone.  Along  with  the  moment  determination  of  M = 11.2  x 10  dyne-cm, 
this  radius  implies  a high  stress  drop  of  about  96  bars.  Evidence  in  the 
amplitude  data  indicates  that  the  polarization  angle  of  shear  waves  is  very 
sensitive  to  lateral  structure. 

A technique  utilizing  theoretical  wave  forms  was  developed  to  determine 
precise  shear  wave  travel  times.  This  technique  was  applied  to  long-period 
World-Wide  Standard  Seismograph  Network  and  Canadian  network  seismograms  of 
five  large  nuclear  explosions  to  obtain  a surface  focus  shear  wave  data  set 
containing  about  100  travel  times  for  distances  greater  than  30°.  Very 
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little  scatter  is  present  in  the  data  from  Novaya  Zemlya  and  the  Nevada  test 
site,  and  so  a reliable  inversion  to  a lower  mantle  velocity  structure  is 
permitted.  This  velocity  model,  based  on  the  59  travel  times  from  Novaya 
Zemlya,  has  significantly  more  structure  than  earlier  models.  The  model, 

SI,  has  proved  to  be  appropriate  for  free  oscillations  as  well  as  for  travel 
times.  This  model  should  be  useful  in  studying  both  lateral  inhomogeneities 
and  the  mineralogical  composition  of  the  earth's  mantle. 

For  the  period  from  15  March  1975  to  30  September  1976  the  principal 
investigators  were  D.  G.  Harkrider  and  1).  V.  Helmberger  under  the  title 
"Body  and  Surface  Wave  Modeling  of  Observed  Seismic  Events."  The  following 
paragraphs  are  a summary  of  the  semi-annual  technical  report.  The  details 
can  be  found  in  sections  III  through  V. 

In  order  to  extend  teleseismic  body  wave  theory  to  more  generalized 
sources  than  a point  shear  dislocation,  it  is  convenient  to  express  the 
source  in  a multipole  form.  An  expansion  of  an  outgoing  linear  elastic  wave 
field  in  spherical  harmonics  provides  an  equivalent  elastic  source  of  quite 
general  character  and  nearly  any  seismic  source  model  can  be  written  in  this 
form.  As  long  as  the  expansion  is  done  in  a linear  homogenious  region, 
potentials  can  be  used  to  continue  the  solution  in  space  and  time  in  the 
region.  In  Section  III  the  equivalent  source  is  then  embedded  in  a stack 
of  plane  elastic  layers  representing  the  near-source  crustal  structure 
and  expressions  are  derived  for  computing  the  steeply  emergent  body  waves 
existing  at  the  base  of  the  model.  These  displacements  can  then  be  combined 
with  transfer  functions  representing  the  effect  of  the  remainder  of  the 
travel  path  to  compute  theoretical  seismograms  for  the  body  waves  recorded 
in  the  far  field.  Using  this  technique  teleseismic  body  waves  have  been 
calculated  by  other  investigators  for  relaxation  sources,  and  for  non  linear 
finite  difference  source  calculations  such  as  for  a three-dimensional  finite 
difference  simulation  of  a strike-slip  earthquake  faulting  and  an  explosion 
in  an  axisymmetric  finite  tunnel. 

One  of  the  interesting  questions  concerning  the  use  of  matrix  methods 
in  the  source  region  is  whether  the  assumption  of  plane  uniform  layers 
implicity  requires  a source  region  which  is  uniformly  layered  outside  of  the 
local  source  region  in  which  teleseismic  body  waves  are  generated  in 
Section  TV,  both  ray  theory  and  propagator  matrix  techniques  were  used  to 
calculate  synthetic  seismograms  for  a point  dislocation  in  a layered  elastic 
medium.  Both  methods  are  equivalent  and  complementary  for  various  model 
calculations.  Ray  theory  has  the  advantage  of  giving  a direct  physical 
interpretation  to  observed  phases.  Matrix  techniques,  while  exact  in  terms 
of  including  all  multiples,  gives  little  insight  into  the  model  but  does 
serve  as  a check  to  the  ray  calculations.  It  can  also  be  used  as  a tool  for 
calculating  the  effect  of  gradients  in  the  eartli  model  by  approximating 
them  with  many  thin  layers. 

Tn  Section  V,  a data  set  of  forty-one  moderate  and  large  earthquakes  has 
been  used  to  derive  scaling  rules  for  kinematic  fault  parameters.  If  effective 
stress  and  static  stress  drop  ayu?equa^^then  fault  rise  time,  r,  and  fault 
area,  S,  are  related  by  t = 16S  *” / ( 7 n g) , where  fi  is  shear  velocity. 


Fault  Length  (parallel  to  strike)  and  width  (parallel  to  dip)  are 
empirically  related  by  L - 2W.  Scatter  for  both  scaling  rules  is  about 
a factor  of  two.  These  scaling  laws  combine  to  give  width  and  rise  time 
in  terms  of  fault  length.  Length  is  then  used  as  the  sole  free  parameter  in 
a Haskell  type  fault  model  to  derive  scaling  laws  relating  seismic  moment  to 
M (20  sec.  surface  wave  magnitude),  M to  S and  m,  (1  sec  body  wave  magnitude) 
to  M . Observed  data  agree  well  with  l:he  predicted  scaling  relation.  The 
"source  spectrum"  depends  on  both  azimuth  and  apparent  velocity  of  the  phase 
or  mode,  so  there  is  a different  "source  spectrum"  for  each  mode,  rather 
than  a single  spectrum  for  all  modes.  Furthermore,  fault  width  (i.e.  the  two 
dimensionality  of  faults)  must  not  be  neglected.  Inclusion  of  width  leads  to 
different  average  source  spectra  for  surface  waves  and  body  waves.  These 
spectra  in  turn  imply  that  m^  and  M reach  maximum  values  regardless  of 
further  increases  in  L and  seismic  moment.  The  m^iM^  relation  from  this 
study  differs  significantly  from  the  Gutenberg-Richter  relation,  because  the 
G-R  equation  was  derived  for  body  waves  with  a predominant  period  of  about 
5 sec  and  thus  does  not  apply  to  modern  1 sec  m,  determinations.  Previous 
investigators  who  assumed  that  the  G-R  relation  was  derived  from  1 sec  data 
were  in  error.  Finally  averaging  reported  rupture  velocities  yields  the 

relation  V = .72?. 
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II.  Abstracts  of  Publications  Not  Previously  Reported  During  this 
Contract  Period 


Langston,  C.  A.  and  R.  Butler,  "Focal  Mechanism  of  the  August  1,  1975, 

Oroville  Earthquake,"  Bull.  Seismo.  Soc.  Am.,  66,  p.  1111-1120. 

Long  period  teleseismic  P and  S waves  from  the  WWSS  and  Canadian  networks 
are  modeled  to  determine  the  focal  parameters  for  the  main  shock  in  the  Oroville 
earthquake  series.  Using  the  techniques  of  P first  motions,  waveform  synthesis, 
and  phase  identification  the  focal  parameters  are  determined  as  follows:  dip  65°; 

rake  -70  ; strike  180  ; depth  5.5  i 1.5  km;  moment  5.7  + 2.0  x 10^  dyne— cm;  and 
a symmetric  triangular  time  function  3 sec.  in  duration.  This  is  a north-south 
striking,  westward  dipping,  normal  fault  with  a small  component  of  left-lateral 
motion.  The  time  separation  between  the  small  foreshock  and  mainshock  appears 
to  be  6.5  sec.  rather  than  8.1  sec.  as  previously  determined. 

Chung,  W.  Y.  and  H.  Kanamori,  "Source  Process  and  Tectonic  Implications  of  the 
Spanish  Deep  Focus  Earthquake  of  March  29,  '954,"  Earth  and  Planet.  Sci.  Letters, 
in  press,  1976. 

The  source  process  of  the  aeep-fi.  cus  Spanish  earthquake  of  March  29,  1954 
(mj,  = 7.1,  h = 630  km)  has  been  studied  by  using  seismograms  recorded  at  tele- 
seismic  distances.  Because  of  its  unusual  location,  this  earthquake  is  considered 
to  be  one  of  the  most  important  earthquakes  that  merit  detailed  studies.  Long- 
period  body  wave  records  reveal  that  the  earthquake  is  a complicated  multiple 
event  whose  waveform  is  quite  different  from  that  of  usual  deep  earthquakes. 

The  total  duration  of  P phases  at  teleseismic  distances  is  as  long  as  40  seconds. 
This  long  duration  may  explain  the  considerable  property  damage  in  Granada  and 
Malaga,  Spain,  which  is  rather  rare  for  deep  earthquakes.  Using  the  azimuthal 
distribution  of  the  differences  between  the  arrival  times  of  the  first,  the 
second  and  later  P phases,  the  hvpocenters  of  the  later  events  are  determined 
with  respect  to  the  first  event.  The  focus  of  the  second  event  is  located  on 
the  vertical  nodal  plane  of  the  first  shock  suggesting  that  this  vertical  plane 
is  the  fault  plane.  This  fault  plane  which  strikes  in  N 2°  E and  dips  89.1°  E 
defines  a nearly  vertical  dip-slip  fault,  the  block  to  the  west  moving  downwards. 
The  time  interval  and  spatial  separation  between  the  first  and  the  second  events 
are  4.3  sec  and  19  km  respectively,  giving  an  apparent  rupture  velocity  of  4.3 
km/sec  which  is  about  74%  of  the  S wave  velocity  at  the  source.  A third  event 
occurred  about  8.8  sec  after  the  first  event  and  about  35.6  km  from  it.  At 
least  6 to  10  events  can  be  identified  during  the  whole  sequence.  The  mechanism 
of  some  of  the  later  events,  however,  seems  to  differ  from  the  first  two  events. 
Synthetic  seismograms  are  generated  by  superposition  of  a number  of  point  sources 
and  are  matched  with  the  observed  signals  to  determine  the  seismic  moment.  The 
seismic  moments  of  the  later  events  are  comparable  to,  or  even  larger  than,  that 
of  the  first.  The  total  seismic  moment  is  determined  to  be  7 x 10-7  dvne-cm 
while  the  moments  of  the  first  and  the  second  shocks  are  2.1  x 10-h  dyne-cm  and 
5.1  x 10-78  dvne-cm  respectively.  The  earthquake  may  respresent  a series  of 
fractures  in  a detached  piece  of  the  lithosphere  which  sank  rapidly  into  the 
deep  mantle  preserving  the  heterogeneity  of  material  property  at  shallow  depths. 
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Source  Parameters  and  Magnitudes 


II 


A data  set  of  41  modertate  and  large  earthquakes  has  been  used  to  derive 
scaling  rules  for  kinematic  fault  parameters.  If  effective  stress  and  static 
stress  drop  are^e^ual,  then  fault  rise  time,  x,  and  fault  area,  S,  are  related 
by  t = 16S  (7tt  6)  , where  6 is  shear  velocity.  Fault  length  (parallel  to  strike) 
and  width  (parallel  to  dip)  are  empirically  related  by  L = 2W.  Scatter  for  both 
scaling  rules  is  abut  a factor  of  two.  These  scaling  laws  combine  to  give  width 
and  rise  time  in  terms  of  fault  length.  Length  is  then  used  as  the  sole  free 
parameter  in  a Haskell  type  fault  model  to  derive  scaling  laws  relating  seismic 
moment  to  Mg  (20-sec  surface-wave  magnitude),  Mg  to  S and  (1-sec  body-wave 
magnitude)  to  Mg.  Observed  data  agree  well  with  the  predicted  scal’’'.g  relation. 
The  "source  spectrum"  depends  on  both  azimuth  and  apparent  velocity  of  the  phase 
or  mode,  so  there  is  a different  "source  spectrum"  for  each  mode,  rather  than  a 
single  spectrum  for  all  modes.  Furthermore,  fault  width  (i.e.  the  two  dimension- 
ality of  faults)  must  not  be  neglected.  Inclusion  of  width  leads  to  different 
average  source  spectra  for  surface  waves  and  body  waves.  These  spectra  in  turn 
imply  that  m^  and  Mg  reach  maximum  values  regardless  of  further  increases  in  L 
and  seismic  moment.  The  mj,  versus  Mg  relation  from  this  study  differs  signifi- 
cantly from  the  Gutenberg-Richter  (G-R)  relation,  because  the  G-R  equation  was 
derived  for  body  waves  with  a predominant  period  of  about  5 sec  and  thus  does 
not  apply  to  modern  1-sec  mg  determinations.  Previous  investigators  who  assumed 
that  the  G-R  relation  was  derived  from  1-sec  data  were  in  error.  Finally, 
averaging  reported  rupture  velocities  yields  the  relation  v^  = 0.726. 


Hart,  R.  S. , D.  L.  Anderson  and  H.  Kanamori,  "Shear  Velocity  and  Density  of  an 
Attenuating  Earth,"  Earth  and  Planet.  Sc  i . Letters,  32^  25-34,  1976. 

The  dispersion  that  must  accompany  absorption  is  taken  into  account  in 
many  recent  body-wave  investigations  but  has  been  largely  ignored  in  surface- 
wave  and  f ree-oscil lat ion  studies.  In  order  to  compare  body-wave  and  f ree- 
oscillation  data  a correction  must  be  made  to  travel  times  or  periods  to  account 
for  absorption-related  physical  dispersion.  The  correction  depends  on  the 
frequency  and  Q of  the  data  and  can  be  as  high  as  1%  which  is  much  larger  than 

the  uncertainty  of  the  raw  data.  Corrected  toroidal  mode  data  is  inverted  to 

obtain  shear  velocity  and  density  versus  depth.  The  average  shear  velocity 
in  the  upper  600  km  is  ^ 2%  greater  than  obtained  from  the  uncorrected  data. 

The  resulting  shear-wave  travel  times  oscillate  about  the  Jef f reys-Bul len  values 
with  an  average  baseline  of  only  +0.5  second.  Thus,  the  discrepancy  between  body- 
wave  and  free-oscillation  studies  is  eliminated. 

Hill,  D.  P.  and  D.  L.  Anderson,  "A  Note  on  the  Earth  Stretching  Approximat ion  for 

Love  Waves,"  Bull.  Seismo.  Soc.  Am.,  in  press,  1976. 

Earth  flattening  transformations  provide  an  efficient  means  for  computing 
Love  wave  dispersion  and  torsional  normal  mode  frequencies  in  radially  hetero- 
geneous, spherically  symmetric  earth  models.  These  transformations  involve 
simple  algebraic  scaling  factors  applied  to  solutions  for  SH  waves  in  a layered 
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half-space.  They  result  In  considerable  computational  savings  over  solutions 
expressed  directly  in  spherical  geometry.  Several  earth  flattening  trans- 
formations for  SH  waves  are  described  in  the  literature  (Anderson  and  Toksoz, 
1973;  Sato,  1968,  Biswas  and  Knopoff,  1970).  Chapman  (1973)  has  examined  the 
general  class  of  power-law  earth  flattening  transformations  and  their  appli- 
cation to  body  wave  problems.  In  this  note  we  clarify  the  earth-stretching 
approximation  used  by  Anderson  and  Toksoz  (1963). 
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ABSTRACT 

The  radiation  field  exterior  to  any  kind  of  volume  source 
in  a homogeneous  medium  can  be  represented  in  terms  of  an  expan- 
sion in  spherical  harmonics.  Such  an  expansion  then  provides  an 
equivalent  elastic  source  representation  of  quite  general 
character  in  that  nearly  any  proposed  seismic  source  model, 
whether  obtained  using  analytical  or  numerical  (finite  differ- 
ence or  finite  element)  methods,  can  be  written  in  this  form. 

The  compatibility  of  this  equivalent  source  with  currently  used 
source  models,  especially  numerical  models  including  detailed 
computations  of  the  nonlinear  processes  at  the  source,  is  dis- 
cussed. The  equivalent  source  is  then  embedded  in  a stack  of 
plane  elastic  layers  representing  the  near-source  crustal  geology 
and  expressions  are  derived  for  computing  the  steeply  emergent 
body  waves  exiting  the  base  of  the  model.  These  displacements 
can  then  be  combined  with  transfer  functions  representing  the 
effect  of  the  remainder  of  the  travel  path  to  compute  theoretical 
seismograms  for  the  body  waves  recorded  in  the  far  field. 


INTRODUCTION 

A fundamental  objective  of  theoretical  seismology  is 
the  development  of  computational  methods  for  accurately 
simulating  the  propagation  of  seismic  waves  through  the  earth. 
The  last  10  - 15  years  have  seen  intense  activity  aimed  at  the 
development  of  techniques  for  computing  elastic  waves  in 
layered  media.  Widespread  use  is  being  made  of  such  methods 
as  generalized  ray  theory  (e.g.,  Gilbert  and  Helmberger,  1972) 
and  the  reflectivity  method  (e.g.,  Fuchs  and  Muller,  1971) 
which  can  compute  the  far  field  body  waves  in  a spherically 
stratified  elastic  earth  model  to  almost  any  desired  precision. 

Concurrent  with  improvements  in  wave  propagation  tech- 
niques, vastly  improved  methods  for  simulating  the  seismic 
source  (earthquake  or  explosion)  have  also  come  into  exten- 
sive use.  These  include  complex  analytical  source  theories 
such  as  that  of  Archambeau  (1968)  as  well  as  the  use  of 
numerical  finite  difference/finite  element  methods.  The 
numerical  methods,  particularly  the  finite  difference  formula- 
tions, are  now  being  applied  in  attempts  to  directly  simulate 
the  nonlinear  processes  at  the  source  (e.g..  Cherry,  et  al . , 
1976b) . 

Much  of  the  data  through  which  we  view  the  seismic 
source  is  recorded  in  the  far  field  and  an  important  test 
of  any  source  model  is  how  well  it  agrees  with  the  far  field 
ground  motion  data.  Therefore,  it  is  necessary  to  have  a 
bridge  between  the  source  calculations  and  the  elastic  wave 
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propagation  techniques  of  theoretical  seismology.  That  is, 
we  need  an  equivalent  elastic  source. 


An  expansion  of  the  outgoing  elastic  wave  field  in 
spherical  harmonics  provides  an  equivalent  elastic  source  of 
quite  general  character  and  nearly  any  seismic  source  model 
can  be  written  in  this  form.  Harkrider  and  Archambeau  (1976) 
derived  the  expressions  for  computing  surface  waves  for  such 
a source  embedded  in  a stack  of  plane  elastic  layers.  In 
this  paper  the  corresponding  theory  for  the  steeply  emergent 
body  waves  exiting  the  base  of  this  plane  layered  model  is 
given.  The  derivation  is  closely  related  to  that  of  Fuchs 
(1966)  and  Hudson  (1969a, b)  who  treated  a similar  problem 
except  that  their  source  was  assumed  to  be  given  in  terms 
of  elementary  point  forces  and  their  derivatives.  Our  re- 
sults are,  of  course,  completely  equivalent  where  the  source 
representations  coincide. 

In  the  following  section  the  form  of  the  equivalent 
elastic  source  and  its  compatibility  with  analytical  and 
numerical  source  calculations  is  discussed.  Subsequent  sections 
contain  the  derivation  of  the  far  field  body  waves  for  such  a 
source  in  a layered  elastic  medium.  Finally,  we  briefly 
discuss  the  computation  of  theoretical  body  wave  seismograms, 
using  this  formulation  for  the  seismic  source  and  crust  in  the 
source  vicinity  together  with  other  methods  for  simulating  the 
rest  of  the  travel  path.  In  a companion  paper  (Bache  and 
Archambeau,  1976) , the  theoretical  seismogram  calculations  are 
discussed  in  greater  detail  and  a number  of  examples  are  given. 
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EQUIVALENT  ELASTIC  SOURCE 


The  radiation  field  exterior  to  any  kind  of  volume 
source  in  a homogeneous  medium  can  be  represented  in  terms  of 
an  expansion  in  spherical  harmonics.  Archambeau  (1968)  seems 
to  have  been  the  first  to  recognize  the  usefulness  of  this 
fundamental  result  and  to  apply  it  to  geophysical  problems. 
The  expansion  in  spherical  harmonics  gives  a compact  equiva- 
lent elastic  source  representation  of  quite  general  character 
and  nearly  any  proposed  seismic  source  model  can  be  cast  in 
this  form.  A brief  description  of  this  source  representation 
and  its  compatibility  with  commonly  used  source  theories  is 
the  subject  of  this  section. 

The  Fourier  transformed  equations  of  motion  in  a homo- 
geneous, isotropic,  linearly  elastic  medium  may  be  written 


--teK’-fe) 


where  u is  particle  displacement  and  and  kg  are  the 

compressional  and  shear  wave  numbers.  The  Cartesian  poten- 
— (4 ) — 

tials  x and  X are  defined  by 

-(4) 

X = 7 • u , 


X = j V x u , 


and  may  be  easily  shown  to  satisfy  the  wave  equation 
V2x(j)  + *-X(j)  = 0,  j = 1,  2,  3,  4, 
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where  k = k = ui/ct  and  k.  = k0  = u>/6  for  i = 1,2,3.  This 

4 U J.  D 

equation  has  as  a solution  the  following  expansion  in  spherical 
eigenfunctions  (e.g.,  Morse  and  Feshbach  (1953)), 

°°  a 

X(j)  (R/W)  = h^2*  (k^R)  [Alm^  ^ cos  m<J> 

£=0  m=0 

+ (ui)  sin  mtp  | P™(cos0)  , (4) 

(2) 

where  the  h„  are  spherical  Hankel  functions  of  the  second 

kind  and  the  P™  are  associated  Legendre  functions.  The 

vector  R has  as  components  the  spherical  coordinates  R,0,<J). 

Equations  (4) , together  with  (1) , provide  an  elastic 

point  source  representation  of  the  (outgoing)  displacement  field. 

The  values  of  the  multipole  coefficients,  A^^  (uj)  , (u)  , 

j = 1,2, 3,4,  prescribe  the  displacement  field  at  all  points  in 

the  homogeneous  medium  where  (1)  applies.  This  point  source 

representation  can  be  viewed  as  a generalized  form  for  a sum 

of  a monopole  or  center  of  dilatation  (£  = 0) , a dipole  or 

couple  (fc  = 1),  a quadrupole  or  double-couple  (i  = 2),  etc. 

For  example,  a center  of  dilatation  is  represented  by  a single 
( 4 ) 

coefficient  A , while  for  a horizontal  double-couple  the 
o o 

nonzero  coefficients  are  -A  ^ = B ^ 2 ^ = A ^ and  B^. 

21  21  22  22 

Description  of  the  character  of  the  elastic  field  gene- 
rated by  seismic  sources  is,  of  course,  a basic  geophysical 
problem.  For  this  paper  it  is  convenient  to  discuss  seismic 
source  descriptions  in  three  categories: 
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1.  Those  obtained  using  finite  difference/finite 
element  numerical  methods. 

2.  Analytical  source  models  of  relaxation  type. 

3.  Dislocation  source  models. 

With  numerical  methods  one  can  attempt  to  directly  in- 
clude complexities  of  the  source  mechanism  in  a deterministic 
computational  scheme.  For  example,  finite  difference  methods 
have  been  extensively  used  to  compute  the  propagating  shock 
wave  due  to  an  underground  nuclear  explosion  (e.g..  Cherry, 
et  al . , 1974)  . In  this  case  the  nonlinear  behavior  of  the 
rock  under  high  stress  loading  determines  the  character  of 
the  seismic  signal.  If  the  source  region  can  be  assumed  to 
be  embedded  in  a medium  in  which  (1)  applies,  an  equivalent 
elastic  source  of  the  form  (4)  can  be  obtained  from  the  out- 
going displacement  field.  This  is  indicated  schematically  in 
Figure  1.  Briefly,  the  procedure  is  to  monitor  the  outgoing 
displacement  field  or,  alternatively,  the  potentials,  x^# 

A 

on  a spherical  surface  of  radius  R.  Using  the  orthogonality 
of  the  spherical  harmonics,  these  potentials  are  related  to 
the  nultipole  coefficients  by 
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Use  of  this  procedure  for  linking  nonlinear  finite  dif- 
ference source  calculations  with  analytical  wave  propagation 
techniques  was  suggested  to  the  first  author  by  Archambeau 
(1973,  personal  communication),  and  has  since  been  implemented 
for  a number  of  complex  explosion  geometries  (Cherry,  et  al . , 
1975,  1976a),  and  for  a three-dimensional  finite  difference 
simulation  of  stick-slip  earthquake  faulting  (Cherry, 
et  al . , 1976b) . The  number  of  terms  required  for  the  expan- 
sion (4)  to  converge  depends  on  the  symmetry  of  the  source 
radiation  at  frequencies  of  interest.  The  most  elementary  ap- 
plication of  the  method  is  for  one-dimensional  (spherically 
symmetric)  explosion  source  calculations.  For  such  problems 
the  elastic  field  is  often  described  by  a reduced  displacement 
potential  defined  by 


U(R,t)  = 


9R 


y (t-R/a)  i 

R J * 


(6) 


Applying  the  Fourier  transform  and  comparing  to  (1)  together 
with  (4),  it  is  easily  derived  that 


A(4)  (oj)  = 
0 0 


-i  k*  y(oi)  , 


(7) 


which  shows  the  equivalence  between  the  reduced  displacement 
potential  and  the  monopole.  For  more  complex  sources  such 
as  an  explosion  in  an  axisymmetric  tunnel  (Cherry,  et.  al . , 
1975)  or  several  explosions  detonated  simultaneously  (Cherry, 
et  al . , 1976a)  quadrupole  and  higher  order  terms  occur  in 
the  expansion.  When  an  earthquake  source  is  computed,  the 
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leading  term  is,  as  expected,  the  quadrupole  (Cherry,  et^  al. , 
1976b)  . 

Three-dimensional  relaxation  models  of  the  seismic 
source  mechanism  have  been  developed  by  Archambeau  (1968,  1972) 
and  Minster  (1973).  These  authors  present  their  results  in 
terms  of  an  expansion  in  spherical  harmonics.  One  form  of 
Archambeau' s model  has  been  used  in  a number  of  studies  of  in- 
duced tectonic  stress  release  by  underground  nuclear  explo- 
sions (Archambeau  and  Sammis,  1970;  Archambeau,  1972;  Bache, 
1976) . In  this  symmetric  problem  only  the  quadrupole  term 
occurs  in  the  expansion.  When  the  Archambeau/Minster  model 
is  used  to  represent  propagating  ruptures  (earthquakes) , many 
terms  may  be  required  for  the  expansion  (4)  to  converge,  de- 
pending on  fault  length,  rupture  velocity  and  frequency. 

The  most  common  approach  to  earthquake  source  theory 
is  to  assume  that  dislocation  theory  is  applicable  (e.g., 
Haskell,  1964;  Savage,  1966).  Such  theories  generally 
represent  the  source  in  terms  of  a double-couple  with  fre- 
quency content  depending  on  the  dislocation  time  history 
assumed.  Harkrider  ( 1976)  has  given  the  expressions  re- 
lating dislocation  source  theories  to  the  expansion  in 
spherical  harmonics  (4).  For  a horizontal  double-couple 
(normal  strike-slip  fault)  Harkrider  (1976),  Equation  (36), 
gives  the  Cartesian  potentials  as 
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0 = 


-i  "-D  (Ui)  k*  sin  20  ^ 

4TTQC02  a 3 


P2(cos0) 

^ »;2><v» . 


0 = 


r = 


-i  H_DWLk.  COM  !Hp!ih(2)(k  , 

47rpai2  B 3 2 6 

?r,  ,,  P1  (cos0)  . 

i kg  sin0  — h(2)  (kRR)  , 

4TTPU)2  6 3 2 0 


(8) 


0 = 


i H-2I2I  k>  cos  2*  h(2)  (k 

4TTPU)2  B 3 2 B 


where  p and  u are  the  density  and  shear  modulus,  and  D(w) 
represents  the  Fourier  transformed  dislocation  time  history. 
The  implied  dimensions  of  these  potentials  are  per  unit  fault 
area.  The  results  can  be  generalized  to  finite  faults  by  in- 
troducing an  integration  over  the  fault  surface.  For  fault 
models  on  which  the  dislocation  history  is  invariant  and  the 
variation  of  phase  between  dislocation  points  and  observer  is 
negligible  over  the  fault  plane,  (8)  multiplied  by  the  fault 
area  are  the  potentials  for  a normal  strike-slip  fault. 

The  Cartesian  potentials  (0,  1JT)  are  related  to  the  dis- 
placements by 

u = V0  + V x 0 . (9) 

Comparing  to  (1) , we  see  that 


(10) 


Using  (10)  and  comparing  (8)  to  (4),  the  multipole 
coefficients  for  the  normal  strike-slip  dislocation  are: 
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b(4> 


(w) 


in  D(W)  k^ 
12tt  (X+2u)  ' 


A(1)  (w)  = 
2 1 


-i  D (to)  kj 
24n 


b(2) 
2 1 


(10) 


A(3)  (to)  = 
2 2 


-A(1)  (to) 
2 1 


(ID 


Expressions  from  which  the  multipole  coefficients  for  a double- 
couple at  arbitrary  orientation  can  be  obtained  are  given  by 
Harkrider  (1976) , Appendix  A. 

Whether  the  seismic  source  is  modeled  by  numerical 
methods  or  analytical  theories,  the  multipole  coefficients 
provide  a computationally  convenient  equivalent  elastic  source. 
Thus  far,  the  discussion  has  been  restricted  to  the  multipolar 
expansion  with  respect  to  a coordinate  system  fixed  with 
respect  to  the  source.  Minster  (1973)  gives  the  transforma- 
tion matrices  by  which  the  multipole  coefficients  in  a standard 
coordinate  system  (e.g.,  fixed  with  respect  to  the  surface  of 
the  earth)  may  be  obtained  from  any  rotated  system.  Using  these 
results,  multipole  coefficients  can  be  computed  in  a convenient 
source  related  system  and  then  rotated  to  a fixed  geographical 
system. 

In  the  following  section  the  equivalent  elastic  source 
in  the  form  discussed  here  will  be  embedded  in  a multilayered 
medium  and  we  will  derive  the  expressions  for  computing  the 
steeply  emergent,  far  field  body  waves. 


III-ll 

THE  EQUIVALENT  ELASTIC  SOURCE  IN  A MULTILAYERED  HALFSPACE 

A number  of  authors  have  investigated  the  elastic  waves 
radiating  from  point  sources  in  a multilayered  medium  over- 
lying  a homogeneous  halfspace.  Most  of  these  studies  have 
relied  on  Thomson-Haskell  matrix  theory  (Thomson,  1950; 

Haskell,  1953)  as  will  our  derivation.  Harkrider  (1964) 
developed  solutions  for  surface  waves  due  to  elementary  point 
forces  at  depth.  Fuchs  (1966)  derived  transfer  functions 
which  include  the  effect  of  the  layered  crustal  model  on  the 
far  field  P waves  from  three  types  of  sources:  a center  of 

dilatation,  a couple  and  a double-couple.  Hudson  (1969a, b) 
extended  Fuchs'  body  wave  results  to  apply  to  quite  general 
sources  of  finite  extent  and  derived  the  analogous  theory  for 
surface  waves.  However,  in  all  these  theories  the  source 
representation  is  in  terms  of  elementary  point  forces  and 
their  derivatives.  The  representation  in  this  form  of  com- 
plicated sources  (including  terms  of  higher  order  than  the 
double-couple)  at  arbitrary  orientation  would  appear  to  be  an 
arduous  task.  Thus,  the  usefulness  of  these  previous  results, 
especially  for  routine  numerical  computations,  seems  to  be 
limited  by  the  inherent  complexity  of  the  source  representa- 
tion . 

The  multipolar  expansion  discussed  in  the  previous 
section  provides  a unique,  compact  and  convenient  numerical 
representation  for  seismic  sources  of  arbitrary  complexity  and 
orientation.  Harkrider  and  Archambeau  (1976)  have  computed  the 
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surface  waves  for  a source  given  in  this  form  embedded  in  a 
multilayered  medium.  This  required  the  formulation  of  the 
displacement  field  in  terms  of  integrals  over  wave  number,  k. 
The  surface  waves  are  then  given  by  the  residue  contribution 
to  these  integrals.  For  the  body  waves  it  is  necessary  to 
evaluate  the  branch  line  contribution  to  similar  integrals  and 
this  is  the  main  result  presented  here.  Derivation  of  the 
k integrals  in  the  appropriate  form  follows  closely  the  deriva- 
tions of  Harkrider  (1964)  and  Harkrider  and  Archambeau  (1976) 
and  the  notation  is,  for  the  most  part,  the  same.  It  should 
be  pointed  out  that  our  results  are  completely  equivalent  to 
those  of  Hudson  (1969a, b) „ The  difference  is,  of  course,  in 
the  source  representation. 

Consider  a semi-infinite  elastic  medium  made  up  of  n 
parallel,  homogeneous,  isotropic  elastic  layers.  Number  the 
layers  from  1 at  the  free  surface  to  n for  the  underlying 
halfspace.  Place  the  origin  of  a cylindrical  coordinate  sys- 
tem (r,<p,z)  at  the  free  surface  and  denote  the  layer  inter- 
faces by  z i = 1,2,  ...,  n-1.  This  geometry  is  depicted  in 
Figure  2.  Let  (q^,  v^,  vT)  be  the  components  of  the  Fourier 
transformed  displacements  in  the  (r,^,z)  directions  in  the  ith 
layer.  Then,  following  Harkrider  (1964),  the  cylindrical 
potentials  ijj\,  are  defined  by 
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1 1 1 — 1 3 


qi  (r , <(> , z) 


9<t.  32^.  . 3fi. 

i , _ i , 1 l 

3r  Sr3z  r 30 


vi (r,4>,z) 


1 . i 32n 

r 34>  r 3z34> 


312. 

l 

3r 


(12) 


wi  (r , $ , z) 


3*.  32^. 

oTz^-  + 7^  + k6  ’ i 

32  3z2  6i  1 


= 1,2, 


n. 


We  will  subsequently  be  interested  in  an  equivalent 

elastic  point  source  at  a depth  z = h.  Let  the  source  layer 

be  denoted  by  a subscript  s;  that  is,  z . < h < z , . It  is 

s s— 1 

necessary  to  express  the  cylindrical  potentials  ($  , Ip  , fi  ) 

s s s 

in  the  source  layer  in  terms  of  the  spherical  potentials,  » 

from  (4).  The  equivalence  is  given  by  Harkrider  and  Archambeau 
(1976)  and  is  as  follows: 


i,  oo 
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where 


X = 

m 


Bm  = 


oo 

1 V''  ni  — n r ,,  -..m  + il  .(4),  , tm+n„mf  ra  \ 

• 77  2^  (1)  isgn(h-2)]  . 


k3 
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°°  /]^r  \ 
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(14) 


■ 77  E (i)m‘n[sgn(h-z)]m*t  A«)  (-l)"*nP,j(^) 


6 £=m 


k0  £=m 

where  j = 1,2,3  in  the  definitions  of  C ^ , D ^ , and 

m r m ' 


m m 


k = u/c,  k^  = u)/y  , 


= (c2/y2-l)  5 , Y = a or  8 , 


(15) 


c = horizontal  phase  velocity. 

The  coefficients  E^,  Fm  in  the  expression  for  ip  are 


to  the 

p(j) 
m ' 

D(*> 

m 
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and,  taking  ci1*  = D*;^  = 0 for  m > l, 
' ^ m m 


2kE  = - C t2]  - , 

m m+1  m-1  m+1  m-1 

2kF  = + C(1]  + F(2]  . d(2]  , 2 < m < l . 

Since  the  location  of  material  boundaries  depends  only 
on  z,  the  dependence  on  r and  4>  will  be  everywhere  the 
same  as  in  the  source  layer.  Therefore,  separate  the  poten- 
tials in  the  layers  as  follows: 


Jm(kr)dk  , 
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f 1 (4>,z) 
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I9  2 y • • « } n • 


The  potentials  satisfy  wave  equations 


for  which  the  general  solutions  may  be  written 

, . , . - ikr  z t \ ' ikr  z 

•{m)(*,z)  - X{m)  e «i  ♦ A[m)  e ai  , 


, . , , -ikra  z ikr  z 

, (m)  ,,  . **  (m)  8-  . ~(m)  8S 

> ■l(<p,z)  = u>-  1 e 1 + u>>  ' e 1 , 


!)«(*, r)  = e'^V  . e-W'  e^V  . 


In  (18)  and  subsequent  equations  subscripts  i,  i = 1,2,  ...,  n, 
denote  quantities  in  the  ith  layer. 
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Then,  define 


S(m) 
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(19) 


We  will  subsequently  be  interested  in  the  values  of  the  poten- 
tial in  the  halfspace  (i  = n) . Applying  the  radiation  condi- 
tion at  infinity,  (18)  and  (19)  give 


^m)(*,z) 
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Now,  combining  (17)  and  (20),  we  have: 


*n(r,$,z) 
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n J (kr)dk  . 
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In  equations  (22)  we  now  have  the  cylindrical  potentials 

<F  , if  , IT  , in  the  nth  layer  (halfspace)  in  terms  of  a sum  of 
n n n 

(ifl)  * (m) 

integrals  of  the  Sommerfeld  type.  The  coefficients  A^  , , 

^ (n\ ) • 

cn  depend  on  azimuth,  <p , as  well  as  wave  number,  k,  and  in- 
clude the  modification  of  the  source  generated  pulse  by  the 
material  discontinuities  in  the  layered  halfspace.  Following 
Harkrider  (1964,  Eqs.  (62)  and  (122)),  these  coefficients  are 
solutions  of  the  mafix  equations: 
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where  the  and  JL  are  given  by  Equations  (61)  and  (124) 

of  Harkrider  (1964).  The  matrix  A is  defined  by  Harkrider 

SI 

and  is  the  layer  product  matrix  which  gives  the  displacement- 

stress  vector  for  P - SV  motion  at  the  source  depth  in  terms 

of  the  displacement-stress  vector  at  the  surface.  Similarly, 

A^  is  the  transfer  matrix  for  the  displacement-stress 
XS1 

vector  associated  with  SH  motion. 

The  source  terms  in  (23)  are  given  by  Harkrider  and 
Archambeau  (1976)  as  follows: 
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The  quantities  rQ,  rg,  P,  u,  Y,  kg  refer  to  the  layer  in  which 
the  source  occurs,  with  p and  y being  the  density  and  shear 
modulus.  The  coefficients  A^,  Bm,  etc.  are  given  by  (14)  with 
the  following  modification.  The  series  are  separated  into  two 
parts ; 


A = A® 
m m 


+ A 


m 


(26) 


with  the  e superscript  denoting  a series  made  up  of  terms 
with  m + l even  and  the  o superscript  denoting  a similar 
series  from  terms  with  m + l odd. 

Equations  (23)  may  be  viewed  as  a set  of  simultaneous 

, . , ^ . . . . , f (m)  ^ (m)  Mm) 

linear  algebraic  equations  and  solved  for  A , u>  , e 

n n n 

in  terms  of  known  quantities.  The  result  is: 
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where 


Equations  (27)  give  A ^ , u>^m\  e ^ in  terms  of  the 
^ n n n 

multipole  coefficients  specifying  the  source  and  the  Haskell- 
Thomson  layer  matrices.  Then  (22)  together  with  (12)  give 
closed  form  solutions  for  the  displacements  in  the  nth  layer 
or  halfspace.  In  the  following  section  we  evaluate  the  inte- 
grals in  (22)  to  extract  the  solutions  for  the  steeply  emergent 
body  waves  of  primary  interest  here. 
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COMPUTATION  OF  BODY  WAVES 

To  compute  the  displacement  potentials  from  (22),  it 
is  necessary  to  evaluate  integrals  of  the  form 
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(k,to) 


-ikr  z 
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e J (kr)  dk  , 
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(29) 


Y = an  or  Bn  , m = 0,  1,  2,  


For  m = 0 the  dilatational  potential,  T , was  evaluated 
at  large  distances  from  the  source  by  Fuchs  (1966)  using  saddle 
point  methods.  Hudson  (1969a, b)  encountered  integrals  very 
similar  to  those  in  (22)  when  solving  for  the  body  waves  due 
to  a point  source  of  general  form  in  a layered  medium.  As 
mentioned  in  the  previous  section,  Hudson's  solution  is 
analogous  to  that  obtained  here,  with  the  difference  being  in 
the  specification  of  the  source. 

Hudson  solved  for  the  far  field  body  waves  given  by 
integrals  like  (29)  by  using  contour  integration  in  the  com- 
plex plane  and  approximating  the  branch  line  contribution  using 
the  saddle  point  method.  The  details  of  this  integration  may 
be  found  in  the  works  of  Fuchs  and  Hudson  and  will  not  be  repro- 
duced here.  It  will  suffice  to  give  the  results  which  are: 
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where  R2  = r2  + z2.  The  geometry  is  shown  in  Figure  2.  The 
saddle  point  approximation  is  valid  as  long  as  R >>  r;  that  is, 
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as  long  as  the  receiver  is  sufficiently  far  from  the  base  of 
the  stack  of  plane  layers.  Alternatively,  since  r/R  = sinB^, 
this  approximation  is  valid  for  steeply  emergent  body  waves. 

With  the  solution  (30),  we  may  now  write  the  far  field 
body  wave  contribution  to  the  potentials  (22) : 
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The  cylindrical  displacements  (q  , v , w ) then  result  from 

n n n 

substituting  (31)  into  (12).  It  is  more  convenient  to  deal 
with  the  displacement  components  in  spherical  coordinates  where 
they  are  identified  as  the  P,  SV  and  SH  components  of  the  pro- 
pagating wave.  In  the  far  field,  these  components  are  given 
by 
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L 


Then  carrying  out  the  indicated  differentiations  and 

_ i 

retaining  only  terns  of  0(R  ),  the  displacements  are: 


(33) 


These  are  the  far  field  body  waves  propagating  into  the  under- 
lying half space  at  specified  horizontal  phase  velocity  c. 
Equations  (33)  , together  with  (27)  , provide  a straightforward 
computational  algorithm  for  computing  these  displacements. 

A summary  of  the  computational  procedure  is  given  below. 

Assume  that  we  have  a seismic  source  specified  in 
terms  of  multipole  coefficients.  We  then  wish  to  compute  the 
P and  S waves  propagating  into  the  halfspace  at  takeoff  angle 
0 , in  which  case  a different  horizontal  wave  speed  is  re- 
quired for  the  P and  S computations.  Alternatively,  we  wish 
to  compute  the  displacements  propagating  at  a fixed  phase 
velocity,  c. 

Having  the  multipole  coefficients,  c and  the  azimuth 
4> , the  following  steps  are  carried  out  for  each  frequency,  u>. 
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1.  Compute  the  source  layer  potentials  $ , 

1*7  , IT  from  (13  - 16)  . 
s s 

2.  Compute  the  displacement-stress  discontinuity 

quantities  <5U  , 6W  , etc.,  from  (24  - 25)  for 
mm 

each  m. 

3.  Compute  the  layer  product  matrices  A , 

R L 

A , J , J . The  latter  two  matrices  are 
LS1 

independent  of  the  source,  while  the  first  two 
depend  only  on  source  depth. 

4.  Compute  the  source  terms  Y , Z , from 

^ m m m 

(29)  . 

5.  Compute  A^m>  , w ^ , e ^ from  (27). 

^ n n n 

6.  Compute  the  displacement  spectra  at  a selected 
distance  R from  (33) . 

Formulation  of  the  solution  in  terms  of  Haskell-Thomson 
matrices  is  convenient  for  various  modifications  that  may  be 
of  interest.  For  example,  Haskell  (1953)  and  Dorman  (1962) 
show  how  to  modify  the  layer  matrices  to  account  for  the 
presence  of  a fluid  layer  and  this  modification  can  easily 


be  carried  through  the  algebra  leading  to  (23)  and  (27) 
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COMPUTATION  OF  THEORETICAL  SEISMOGRAMS 

Our  objective  is  to  develop  improved  methods  for  com- 
puting theoretical  seismograms  at  large  (say  greater  than 
1500  km)  distances.  Equations  (33)  give  the  P,  SV,  SH 
waves  emanating  into  a homogeneous  half space  underlying 
a stack  of  plane  layers  in  which  a seismic  source  is  em- 
bedded. If  this  plane  layered  model  is  used  to  represent 
the  crust  in  the  source  vicinity,  (33)  represents  the 
waves  propagated  into  the  upper  mantle.  The  major  restriction 
on  the  use  of  (33)  is  to  situations  for  which  the  saddle 
point  approximation  is  valid;  that  is,  to  small  angles  of 
emergence . 

To  compute  theoretical  seismograms  we  need  to  combine 
(33)  with  other  methods  for  computing  the  remainder  of  the 
travel  path  to  the  receiver.  There  are  a number  of  computa- 
tional schemes  that  could  be  chosen  but  we  have  found  it  con- 
venient to  use  generalized  ray  theory  (Wiggins  and  Helnberger, 
1974)  for  the  upper  mantle  and  Haskell-Thomson  matrices 
(Haskell,  1962)  for  the  crust  in  the  vicinity  of  the  receiver. 
Thus  we  break  the  travel  path  into  three  pieces: 

1.  The  crust  at  the  source  down  to  a depth  D. 

2.  The  crust  at  the  receiver  to  the  same  depth  D. 

3.  A laterally  homogeneous  upper  mantle  extending 
from  D to  depths  greater  than  the  deepest  turning 
point  of  interest. 
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The  travel  path  segments  are  linked  by  requiring  that 
the  velocities  at  D be  the  same  in  all  three  structures. 

This  scheme  gives  a great  deal  of  computational  flexibility 
in  that  portions  of  the  travel  path  can  be  varied  without 
repeating  the  entire  calculation.  The  actual  implementation 
is  described  in  a companion  paper  (Bache  and  Archambeau, 

1976)  where  a number  of  examples  are  presented. 

Methods  similar  to  those  outlined  here  have  been  used 
by  a number  of  investigators.  Most  closely  related  is  the 
work  by  Douglas  and  colleagues  (e.g.,  Douglas,  et  al . , 1972, 
1974;  Cullen  and  Douglas,  1975)  which  uses  the  method  of 
Hudson  (1969a, b)  for  embedding  a seismic  source  in  a plane 
layered  model  of  the  source  crustal  structure.  Also,  these 
authors  represent  the  geometric  spreading  effect  of  the  upper 
mantle  by  a constant  which  is  a function  only  of  epicentral 
distance  (Carpenter,  1966)  rather  than  the  detailed  general- 
ized ray  theory  method  we  prefer.  We  note  that  beyond  the 
triplications  the  upper  mantle  effect  does  essentially  re- 
duce to  a distance  dependent  constant,  though  its  value  de- 
pends on  the  earth  model  used.  Other  authors  who  have  used 
similar  techniques  include  Kogeus  (1968)  and  Hasegawa  (1972, 
1973)  who  used  Fuchs  (1966)  formulation  for  the  source  crustal 
transfer  function  and  Julian  and  Anderson's  (1968)  geometric 
spreading  factor. 

For  simple  point  sources  (center  of  dilatation,  couple 
or  double-couple)  generalized  ray  theory  can  be  directly  ap- 
plied to  compute  theoretical  seismograms  (e.g.,  Wiggins  and 


Ill -28 


Helmberger,  1974;  Muller,  1971).  Another  successful  theo- 
retical seismogram  computing  technique  is  the  reflectivity 
method  (Fuchs,  1975;  Fuchs  and  Muller,  1971;  Kind  and  Muller, 
1975)  though  it  too  is  restricted  to  elementary  point  source 
representations . 

The  major  difference  between  the  method  presented  in 
this  paper  and  the  others  mentioned  above  is  in  the  source 
representation.  Our  method  is  essentially  independent  of 
the  complexity  of  the  source  as  long  as  the  outgoing  elastic 
wave  field  can  be  expanded  in  spherical  harmonics. 
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FIGURE  CAPTIONS 

1.  Schematic  display  of  the  determination  of  an 
equivalent  elastic  representation  for  an  arbitrary 
volume  source. 

2.  The  geometry  and  coordinate  system  for  a source  at 
depth  h in  a multilayered  halfspace. 
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Abstract 

Simultaneous  modelling  of  source  parameters  and  local  layered  earth 
structure  for  the  April  29,  1965,  Puget  Sound  earthquake  was  done  using  both 
ray  and  layer  matrix  formulations  for  point  dislocations  imbedded  in 
layered  media.  The  source  parameters  obtained  are:  dip  70°  to  the  east, 

9 6 

strike  344°,  rake  -75°,  63  km  depth,  average  moment  of  1.4  + 0.6  x 10  dyne-cm, 
and  a triangular  time  function  with  a rise  time  of  0.5  sec  and  fall-off 
of  2.5  sec.  An  upper  mantle  and  crustal  model  for  southern  Puget 
Sound  was  determined  from  inferred  reflections  from  interfaces  above  the 
source.  The  main  features  of  the  model  include  a distinct  15  km  thick  low 
velocity  zone  with  a 2.5  km/sec  P wave  velocity  contrast  lower  boundary 
situated  at  approximately  56  km  depth.  Ray  calculations  which  allow  for 
sources  in  dipping  structure  indicate  that  the  inferred  high  contrast  value 
can  trade  off  significantly  with  interface  dip  provided  the  structure  dips 
eastward.  The  crustal  model  is  less  than  15  km  thick  with  a substantial 
sediment  section  near  the  surface.  A stacking  technique  using  the 
instantaneous  amplitude  of  the  analytic  signal  is  developed  for  interpreting 
short  period  teleseismic  observations.  The  inferred  reflection  from  the  base 
of  the  low  velocity  zone  is  recovered  from  short  period  P and  S waves.  An 
apparent  attenuation  is  also  observed  for  pP  from  comparisons  between  the 
short  and  long  period  data  sets.  This  correlates  with  the  local  surface 
structure  of  Puget  Sound  and  yields  an  effective  Q of  approximately  65 
for  the  crust  and  upper  mantle. 
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Introduction 

A previous  naoer  dealt  with  the  problem  of  trying  to  deduce 
source  parameters  from  an  extremely  shallow  earthquake  . the  Koyna  event  of 
December,  1967  (Langston,  1976a).  The  interference  of  the  direct  waves  and  surface 
reflections  was  very  severe  due  to  the  depth  of  the  Koyna  hypocenter.  S>me  speculations 
were  made  to  the  origin  of  some  of  the  reverberations  after  the 
reflections  but  these  were  hampered  by  both  ignorance  of  the  Koyna 
crustal  structure  and  the  interference  with  the  major  phases.  If, 
however,  an  earthquake  is  deep  enough  so  that  the  surface  reflections 
and  the  direct  wave  are  well  separated,  perhaps  layer  interfaces 
between  the  hypocenter  and  free  surface  can  be  resolved  by  intermediate 
reflected  arrivals.  It  is  exactly  this  supposition  which  will  be 
used  to  explain  the  shape  of  long  period  P waves  from  the  Puget  Sound 
earthquake.  A crust  model  for  southern  Puget  Sound  will  be  deduced 
by  identifying  arrivals  between  P and  pP  as  underside  reflections  from 
crustal  layers.  A major  compressional  and  shear  wave  low  velocity 
zone  (LVZ)  in  the  uppermost  mantle  will  be  proposed  by  modelling  the 
time,  amplitude,  and  polarity  of  these  reflected  arrivals. 
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Simultaneously,  care  will  be  taken  to  model  the  source  parameters  as  well. 
Geologic  and  Tectonic  Sot t ing 

Puget  Sound  is  located  in  northwestern  Washington  and  is  part  of  a 
major  north-south  geologic  and  physiographic  province,  the  Puget-Wil lamette 
depression  (Figure  1)  (Snavely  and  Wagner,  1963).  The  geologic  development 
of  this  area  has  been  controlled  by  Tertiary  sedimentation  and  volcanism 
in  a north-south  trending  eugeosyncline  shown  in  Figure  1 (Snavely  and 
Wagner,  1963;  Snavely  _et  a_l.  , 1968).  The  Puget-Willamette  depression  may 
be  fault  controlled  with  major  faults  occurring  on  the  eastern  margin  of 
Puget  Sound  and  the  western  side  of  the  Willamette  Valley  (Algermissen  and 
Harding,  1965;  Bromery  and  Snavely,  1964). 

Of  seismic  refraction  profiles  which  have  been  made  on  land,  perhaps  the 
most  significant  is  the  one  by  Berg  et_  a_l . (1966).  They  obtain  a very  thin 
(16  km)  crust  and  a Pn  velocity  of  8.0  km/sec  for  the  Coast  Ranges  of  Oregon 
and  southern  Washington.  fatal  and  Tuve  (1955),  who  did  some  of  the  first  crustal 
refraction  studies  in  the  Pacific  northwest,  reached  a similar  conclusion 
for  the  southern  Puget  Sound  area  obtaining  a crustal  thickness  of  ibout  19  km 
and  Pn  velocity  of  8.0  km/sec.  Because  of  the  conspicuous  lack  of  pre-Tertiary 
granitic  basement  within  the  margins  of  the  Tertiary  eugeosyncline  (Figure  1) 
and  from  the  velocity  and  thickness  values  obtained  in  this  refraction 
profile  some  authors  have  suggested  that  the  area  is  a large  cmhuvment  of 
oceanic  crust  in  the  North  American  continent  (Hamilton  and  Myers,  196b; 


Dickinson,  1970). 


IV  - 3. 


Through  detailed  gravity  and  magnetic  modelling  Danes,  et  al.  (1963)  deduced 
some  of  the  major  structures  of  southern  Puget  Sound.  Their  principal 

results  include  the  discovery  of  a large  northwest  trending  igneous 
horst  structure  with  flanking  deep  sedimentary  basins.  They 
infer  sediment  thicknesses  of  4 km  and  10  km  for  the  southern  Tacoma 
basin  and  northern  Seattle  basin,  respectively.  The  large  observed 
gravity  gradients  suggest  major  faults  bounding  this  horst  with 
inferred  throws  on  the  order  of  several  kilometers.  It  may  be 
significant  that  the  April  1965  event  occurred  under  the  eastern 
edge  of  this  structure.  Danes  (1969)  has  also  suggested  that  Puget 
Sound  is  a great  isostatic  depression  with  an  isostatic  anomaly 
greater  than  -100  mgal. 

The  theory  of  plate  tectonics  forms  the  basis  for  current  thinking 
on  the  tectonic  framework  of  the  Pacific  Northwest.  Interpretation 
of  the  magnetic  anomalies  of  Raff  and  Mason  (1961)  by  Vine  and  Wilson 
(1965),  Wilson  (1965),  and  Vine  (1966)  and  other  subsequent  work  by 
Tobin  and  Sykes  (1968),  Atwater  (1970),  Dickinson  (1970),  Silver 
(1971a,  1971b,  1972)  , Crosson  (1972),  and  Chandra  (1974),  among  others, 
have  led  to  a hypothesis  relating  the  geophysical  and  geologic  data 
into  one  plate  tectonic  scheme.  The  model  states  that  subdue tion  of 


oceanic  crust  and  upper  mantle  is  taking  place  north  of  Cape 
Mendocino  to  Vancouver  Island.  A small  offshore  ridge  system,  the 
Gorda-Juan  de  Fuca  rise,  is  still  producing  oceanic  lithosphere 
which  is  the  remnant  of  the  previously  more  extensive  Farallon  plate 
(Atwater,  1970).  Spreading  at  a half-rate  of  about  2.5  cm/year, 
this  small  plate  may  still  be  underthrusting  the  continent  as  the 
association  of  andesitic  vulcanism  in  the  Cascades  seems  to  require 
(Dickenson,  1970)  and  as  inferred  from  offshore  geologic  evidence 
(Silver,  1971a, 1972).  Against  the  normal  assumptions  of  plate 
tectonics,  this  small  plate  does  not  appear  to  be  rigid,  but  seems  to 
be  experiencing  internal  deformation  (Silver ,1971b ; Crosson,  1972). 

Perhaps  the  most  serious  problem  with  this  scheme  is  the 
conspicuous  absence  of  a seismic  Benioff  zone.  The  general  seismicity 
level  of  the  area  is  low  other  than  at  the  offshore  fracture  zones 
(Tobin  and  Sykes,  1968;  Chandra,  1974).  Puget  Sound  has  a moderate 
background  of  diffuse  seismicity,  however,  and  a few  events  have 
occurred  at  depths  of  up  to  70  km  (Crosson,  1972).  The  magnitude 
6.5  Puget  Sound  earthquake  of  April  1965  was  located  at  60  kra  depth, 
see  Figure  1 (Algermissen  and  Harding,  1965),  and  the  somewhat 
larger  April  1949  event,  at  70  km  (Nuttli,  1 0 5 2 ) . Although  these 
events  are  not  extremely  deep  for  other  island  and  continental  arcs, 
they  are  very  unusual  compared  to  all  other  continental  U.  S. 
earthquakes.  Because  of  this,  they  have  been  cited  ns  evidence  for 


rv- 


a remnant  of  a subducting  plate,  or  perhaps,  a very  slowly  subducting 
plate  under  Washington  (Isacks  and  Molnar,  1971;  Crosson,  1972). 

Other  geophysical  evidence  supporting  this  plate  model  comes  from 
travel  time  anomalies  determined  from  P residuals  for  the  1965  event 
(McKenzie  and  Julian,  1971)  and  from  an  array  processing  study  using 
teleseismic  P arrivals  at  Puget  Sound  (Lin,  1974).  Both  of  these 
studies  reach  the  similar  conclusion  that  a high  wave  velocity  plate 
dips  to  the  east  with  an  angle  of  about  50  degrees. 

The  unusual  depth  of  the  April  1965  event  and  the  interesting 
geologic  and  geophysical  problems  this  area  presents  motivates  the 
detailed  waveform  study  that  this  paper  reports.  The  original 
purpose  of  this  study  was  to  find  a detailed  source  model  for  the 
earthquake,  but  as  it  turned  out,  much  more  interesting  information 
could  be  found  by  modelling  the  local  source  structure  as  well. 

The  Puget  Sound  Earthquake 

On  April  29,  1965,  at  15:28:44  GMT  a magnitude  6.5  earthquake 
shook  the  environs  of  southern  Puget  Sound  causing  moderate  damage 
in  Seattle  and  Tacoma.  The  location  of  the 

epicenter  was  midway  between  these  two  cities  and  the  hypocontral 
depth  estimated  to  be  at  60  km  (Algermissen  and  Harding,  1965). 
Consistent  focal  mechanisms  for  the  event  done  by  several  authors 
show  predominantly  normal  dip-slip  movement  on  a 70°  eastward 
dipping  plane  striking  approximately  15°  west  of  north 
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(Algermissen  and  Harding,  1965;  Isacks  and  Molnar,  1971;  Chandra, 
1974).  The  auxiliary  plane  is  only  poorly  constrained  because  of 
sparse  local  station  coverage,  a common  occurrence  for  this  type  of 
orientation.  Because  of  excellent  teleseismic  coverage,  however, 
the  first  nodal  plane  is  extremely  well  determined  and  serves  as  a 
very  useful  constraint  in  the  waveform  modelling. 

Data  and  Data  Processing 

The  gathering  and  processing  of  long  period  P and  SH  waveforms 
were  done  as  described  in  Langston  (1976a) • Table  1 
lists  the  WWSSN  stations  utilized  for  this  study.  Unfortunately, 
there  was  only  one  station  in  which  the  horizontal  components  were 
naturally  rotated  with  respect  to  the  ray  direction.  As  a result, 
only  a few  SH  waveforms  could  be  salvaged  for  waveshape  analysis 
and  even  these  may  be  contaminated  by  the  rotation  process. 

Most  of  the  stations  in  Table  1 were  equipped  with  the  longer 
period  30-100  instrument  instead  of  the  15-100  instrument  used  in 


the  Koyna  study.  A few  stations,  ANT,  QUI,  SJG,  and  HF.C  had  just 
been  changed  to  the  standard  15-100  instrument,  however,  so  these 
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were  equalized  to  be  consistent  with  the  rest  of  the  data  set.  An 
operator  was  constructed  and  convolved  with  these  data  to  effectively 
make  them  30-100  observations.  The  30-100  instrument  response  was 
calculated  using  Hagiwara's  (1958)  formulation. 

Data  Inversion  and  Interpretation 

As  a starting  point,  it  would  be  most  convenient  to  present  the 
final  inferred  earth  and  source  models  and  P waveform  fits.  The 
complex  interactions  of  the  earth  and  source  models  will  then  be 
examined  point  by  point  and  the  reasoning  behind  each  effect  presented 
in  a coherent  manner. 

Figure  2 displays  the  results  of  trial  and  error  waveform 
modelling  to  find  a source  and  earth  model  most  consistent  with  the 
long  period  data.  A standard  first  motion  plot  (bottom  hemisphere) 
is  given  in  the  center  of  the  diagram  with  the  P nodal  planes 
inferred  from  this  study.  In  this  determination,  pP  as  well  as  P 
first  motions  were  incorporated.  Lines  are  drawn  from  each  observed- 
synthetic  waveform  pair  to  the  appropriate  spot  on  the  focal  sphere 
which  represents  the  downgoing  P ray  at  that  station.  For  each 
seismogram  pair  the  observed  is  plotted  above  the  synthetic.  The 
source  orientation  parameters  are  given  in  the  corner  of  Figure  ? 

In  this  model  the  source  is  assumed  to  be  a single  point  dislocation 
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with  a triangular  time  function  characterized  by  a rise  and  fall-off  time 
St^  and  <5t2  equal  to  0.5  and  2.5  sec,  respectively.  The  final  plane- 
layered earth  model  is  presented  in  Figure  3 and  is  under  the  heading  of 
'PS-9'  in  Table  2.  The  source  is  situated  at  a depth  of  63  km  for  this  model. 

Both  ray  theory  (Langston  and  Helraberger,  1975)  and  propagator  matrix 
techniques  (Haskell,  1953;  Harkrider,  1964;  Fuchs,  1966)  were  used  to  calculate 
synthetic  seismograms  for  a point  dislocation  in  a layered  elastic  medium. 

Both  methods  are  equivalent  and  complementary  for  various  model  calculations. 

Ray  theory  has  the  advantage  of  giving  a direct  physical  interpretation  to 
observed  phases.  Matrix  techniques,  while  exact  in  terms  of  including  all  multiples, 
gives  little  insight  into  the  model  but  does  serve  as  a check  to  the  ray  calculations. 
It  can  also  be  used  as  a tool  for  calculating  the  effect  of  gradients  in  the  earth 
model  by  approximating  them  with  many  thin  layers. 

Two  striking  features  are  apparent  in  Figures  2 and  3.  The 
first  is  the  quality  of  the  waveform  fits  over  the  entire  azimuth 
range.  The  observed  waveforms  have  a tremendous  variation  in  shape 
as  a function  of  azimuth  and  take-off  angle  which  the  model  approximates 
quite  nicely  for  nearly  all  stations.  Secondly,  the  earth  model 
presented  in  Figure  3 is  quite  unusual.  The  major  features  of  PS-9 
include  a very  distinct  and  large  low-velocity  zone  between  41  and 
56  km  depth.  This  structure  is  sandwiched  between  what  appears  to  bo 
mantle  velocity  material.  The  crustal  section  at  the  top  is  very  thin, 
less  than  15  km,  and  has  very  low  velocity  materials  near  the  free 
surface.  This  model  was  inferred  entirely  from  the  long  period  P waves 
and  will  be  discussed  by  closely  examining  which  characteristics  of 


the  P waves  control  its 


various  detni 
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Let  us  first  look  at  what  effects  a simple  earth  model,  a layer 
over  a halfspace,  gives  for  the  long  period  P response.  Figure  4 
compares  the  simple  one  layer  model  (Table  2 ) with  PS-9  for 
representative  P waveforms.  The  major  phases  P,  pP,  and  sP  are 
quickly  apparent  in  these  waveform  comparisons  although  there  are 
significant  differences  in  the  interval  between  the  direct  P wave 
and  pP.  The  small  arrival  several  seconds  before  pP  in  the  one  layer 
model  is  the  P reflection  from  the  bottom  of  the  layer  except  at  GUA 
where  it  is  an  b-P  conversion  at  reflection.  The  first  18  sec  of  the 
P-wave  for  the  simple  model  is  just  the  convolution  of  the  source  time 
function  with  the  Q operator  and  the  instrument 

response  since  there  are  no  distortions  due  to  near  source  structure. 

For  COP  and  stations  like  it  in  the  northeast  (see  Figure  2) 
the  first  18  sec  of  the  observed  P wave  show  at  least  two  arrivals 
after  the  dilatational  direct  P.  First  there  is  an  equal  sized 
compression,  arrival  A (Figure  4)  approximately  two  seconds  after 
P with  another  dilatational  arrival,  arrival  B,  2 seconds  after 
that.  F.xamination  of  the  one  layer  model  of  COP  demonstrates  that 
the  overswing  of  the  instrument  response  is  not  a factor  here.  At  PTO 
the  direct  wave  and  arrival  A have  the  same  sign  and  add  constructively 
because  of  the  change  in  the  P radiation  pattern.  This  shows  up  as  an 
increase  in  the  P amplitude,  relative  to  pP,  at  PTO  and  at  other  similar 
stations.  Arrival  B is  plainly  the  same  polarity  as  at  COP.  The 
interpretation  of  these  arrivals  is  based  on  the  assumption  that 
they  represent  underside  reflections  from  discontinuities  above  the 
source  and  not  from  source  complications.  It  is  entirely  reasonable 
to  assume  that  there  are  major  discontinuities  between  a 60  km  deep 
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source  and  the  free  surface, but  it  is  obviously  extremely  hard  to 
prove, unequi vocally,  that  small  arrivals  are  from  such  reflections 
and  not  source  effects.  The  seismograms  of  Figure  2 suggest  that 
these  arrivals  are  associated  with  pP  rather  than  the  direct  P wave 
by  the  observation  that  when  direct  P changes  polarity  the  A and  B 
arrivals  remain  constant. 

Assuming  that  these  arrivals  came  from  reflectors  above  the 
source,  what  can  be  said  about  their  properties?  Since  the  strength 
of  up-going  P determines  the  amplitude  of  the  reflection  as  well  as 
the  material  contrasts,  an  approximate  determination  of  the  velocity 
contrast  can  immediately  be  made.  At  the  northeastern  P nodal 
stations  upgoing  SV  radiation  is  at  minimum  because  of  fault 
orientation  so  that  only  P interactions  can  be  considered.  A 
cursory  examination  of  the  waveforms  gives  a minimum  amplitude 
estimate  for  reflection  A of  about  0.15.  This  is  clearly  an  under- 
estimate since  interference  with  the  direct  wave  knocks  it  down 
considerably.  There  must  be  interference  since  the  width  of  the 
first  swing  changes  as  a function  of  P amplitude  and  azimuth.  Compare 
STU  and  T01,,  Figure  2,  for  example.  The  type  of  contrast  can 
immediately  be  deduced  because  of  the  known  polarity  of  pP.  Since 
upgoing  P is  dilatational  and  the  reflection  is  compi essional , the 
reflection  coefficient  must  be  negative,  which  implies  a higher 
to  lower  velocity  contrast.  Simultaneously  modelling  the  time  function, 
relative  timing,  and  relative  amplitudes  of  direct  P and  phase  A for 


the  northeastern  and  southeastern  stations  yields  the  high  contrast 
of  8.0  to  5.5  km/sec  for  the  lowermost  boundary  of  PS-9,  Figure  3. 

Continuing  this  line  of  reasoning  one  step  further  to  phase  B 
gives  some  remarkable  results.  Using  the  same  process  for  finding  the 
sign  of  the  reflection  coefficient,  the  polarity  of  B suggests 
that  it  comes  from  a lower  to  higher  velocity  contrast.  Phase  B, 
therefore,  delimits  the  top  of  a low  velocity  zone.  Figure  5 
demonstrates  this  explicitly  for  the  station  KEV.  Ray  number  one  is 
normalized  to  unity  and  all -other  ray  amplitudes  referenced  to  it. 

The  top  of  the  LVZ  was  modelled  with  two  in~te~r i a ce s in  order  to 
increase  the  width  of  the  reflected  pulse.  This  particular  model 
explains  the  azimuthal  variation  of  wave  shapes  for  the  first  eighteen 
seconds  remarkably  well  (Figure  2). 

The  uncertainty  in  orientation  (see  Figure  2)  is  not  a major 
factor  in  the  modelling.  The  T first  motion  data  constrain  the 
north-south  nodal  plane  to  within  a degree.  Since  all  of  the  stations 
are  near  the  center  of  the  plot  and  most  are  close  to  the  nodal 
pi ane, relat ive ly  large  variations  in  the  rake  (-90  ±20°)  do  little 

to  affect  the  relative  amplitudes  of  P and  pP.  Upgoing  SV  is  more 
sensitive  so  there  is  some  basis  for  assigning  the  particular  value 


used,  although  it  is  a small  effect  and  will  be  discussed  later. 

The  velocity  assigned  to  the  source  layer  is  a v.ilue  typically  found 
for  these  depths  in  upper  mantle  studies.  Variations  in  this  velocity, 
of  up  to  * 0.5  km/sec,  don't  significantly  affect  the  results  since 
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we  are  looking  at  relative  amplitudes  and  velocity  contrasts  only. 

This  does  point  out  that  the  absolute  values  of  velocity  for  any 
part  of  the  inferred  model  are  somewhat  ambiguous. 

The  phase  pP  and  sP  control  and  constrain  the  top  of  model  PS-9. 

For  the  source  function  inferred  from  the  direct  P wave,  it  is 
evident  that  pP  in  a simple  one-layer  model  (Figure  4)  starts  too 
sharply.  The  small  arrival  before  pP  suggests  that  a number  of  small 
reflections  in  the  upper  crust  could  give  the  desired  effect  of 
producing  a smooth  ramp  before  the  main  pP  peak.  Further  evidence 
for  this  type  of  model  occurs  in  the  western  stations  HNR,  GUA,  and 
ANP . The  predicted  sP  phase  for  the  one  layer  model  is  much  too 
large  and  is  not  affected  by  small  changes  in  the  radiation  pattern. 

The  easiest  way  to  cut  down  this  amplitude  is  to  lower  the  reflection 
coefficient  by  using  low  velocities  near  the  top  of  the  model.  This 
implies  many  contrasts  also,  since  there  must  be  some  kind  of 
transition  from  mantle  to  sediment  velocities.  This  is  the  basis  for 
modelling  the  upper  crust  in  model  PS-9.  It  is  interesting  to  note 
that  the  crustal  section  had  to  be  kept  thin  since  thicker  crustal  models 
caused  spurious  arrivals  from  the  Moho  before  the  arrival  time 
of  pP.  The  absolute  velocities  for  the  top  are  not  to  veil  constrained. 
These  values  were  determined  by  using  the  local  refraction  results 
of  Berg  et  al.  (1966)  and  Tatel  and  Tuve  (1955).  The 
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lowest  velocities  are  appropriate  for  Tertiary  clastic  rocks  as 
reported  in  Press  (1966).  The  mantle  velocity  above  the  LVZ  of  7.8 
km/sec  was  based  on  the  data  of  Dehlinger  et  JalL.  (1965)  and  McCollum 
and  Crosson  (1975). 

The  S wave  velocity  is  one  of  the  least  constrained  parameters 
in  the  model.  The  average  S time  is  only  constrained  by  sP,  a 
phase  which  is  distorted  and  reduced  in  size  by  the  model.  However, 
there  is  some  evidence  that  large  S wave  contrasts  exist  in  the 
LVZ.  At  the  western  azimuth  stations  S-P  reflections  were  needed 
to  reduce  the  backswing  of  the  direct  P (see  GUA,  Figure  4, for  a 
comparison)  and  theoretical  arrivals  after  sP  were  only  obtained 
after  increasing  the  shear  velocity  contrast  at  the  boundaries  of 
the  LVZ.  These  arrivals  are  shown  as  'C'  and  ' D*  in  Figure  4 and 
are  crustal  multiples  with  substantial  S to  P conversions.  Figure  2 
shows  that  these  multiples  can  help  explain  the  arrivals  after  sP 
at  these  stations. 

These  multiples  bring  up  interesting  questions  concerning  the 
arrivals  after  pP  at  the  northeastern  and  eastern  stations.  PS-9 
predicts  a few  crustal  multiples  after  pP  but  nothing  like  those 
in  some  of  the  observations.  Figure  6 shows  the  long  period  P 
wave  from  the  station  KEV  with  a synthetic  produced  from  the  PS-9 
model  above  and  one  made  from  a preliminary  model,  PS-1  (Table  2), 
below.  The  PS-1  model  did  not  predict  the  LVZ  interference  effects 
very  well  nor  the  shape  of  pP  for  most  observations.  Because  of  the 
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low  shear  velocities  and  higher  contrasts  this  model  has, compared 
to  PS-9,  crustal  reverberations  after  pP  are  quite  pronounced  and 
fit  the  observations  quite  nicely  for  this  azimuth.  However, 
since  it  did  not  fit  the  front  part  of  the  record  very  well,  in 
general,  it  was  not  used.  This  exercise  demonstrates  that  the 
arrivals  after  pP  could  conceivably  be  explained  by  large  velocity 
contrasts  although  lateral  inhomogeneity  would  probably  be  needed 
to  match  them.  These  crustal  and  mantle  reverberations  sample 
larger  portions  of  the  model  and  at  distances  further  away  from  the 
epicenter  as  relative  arrival  times  increase.  Lateral  changes  over 
a scale  length  of  50  km,  not  inconceivable  for  the  region,  coupled 
with  the  substantial  depth  of  the  source  could  be  responsible  for 

such  effects.  It  must  be  mentioned  that  an  added  ambiguity  inherent 

in  this  modelling  is  the  lack  of  receiver  characteristics.  For  the 
same  reasons  cited  in  the  Koyna  study  (Largston,  197ba)  no  receiver 
responses  were  evaluated.  Presumably,  the  effect  for  the  vertical 
P wave  is  small  but  could  be  on  the  same  order  as  the  small  arrivals 

behind  pP.  As  such,  these  unknowns  have  to  be  considered  a source 

of  error  in  the  study. 

The  evidence  for  constraining  the  rake  angle  of  the  north-south 
nodal  plane  comes  from  the  relative  amplitude  of  the  reflection  from 
the  bottom  of  the  LVZ.  The  observations  of  Figure- 3. 2 suggest  that 
the  reflection  is  more  pronovinced  for  the  northern  stations 
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(REV,  NUR,  etc.)  compared  to  the  southern  stations  (QUT,  NNA,  LPB, 

ANT).  Rather  than  invoke  extreme  lateral  heterogeneity  for  this 
interface,  a simple  explanation  can  be  made  for  the  effect  using 
fault  orientation.  If  the  fault  was  pure  dip-slip  (A  = - 90°)  the 
northern  and  southern  stations  should  be  identical  since  the  radiation 
pattern  would  be  symmetric  about  a line  perpendicular  to  the  fault 
plane.  Varying  the  rake  angle  strongly  affects  only  the  upgoing  SV 
radiation  since  it  is  near  a node  for  eastern  azimuths  and  this  type 
of  orientation.  It  was  found  that  the  interference  of  S-P  conversions 
with  the  P reflections  decreased  the  amplitude  of  the  LVZ  phases 
where  upgoing  SV  was  comparatively  large.  Using  this  effect  the  rake 
angle  was  deduced  to  be  approximately  -75°  rather  than  pure  dip-slip. 

Because  the  extremely  high  velocity  contrasts  found  for  the  LVZ  of 
PS-9,  it  might  prove  interesting  to  perform  a parameter  study  which  includes 
one  more  level  of  structure  complication  by  including  layer  dip.  Figure  7 
is  the  result  of  such  a calculation.  The  bottom  interfaces  of  PS-9  which 
include  the  LVZ  were  simply  tilted  towards  the  east  preserving  all  laver 
thicknesses  and  velocities.  A ray  theoretical  approach  was  used  to  calculate 
the  response  for  this  three  dimensional  model  and  is  described  elsewhere 
(Langston,  1976c).  Since  all  interfaces  are  parallel,  the 

method, in  this  case,  simply  consists  of  assigning  each  ray  an  amplitude 
corresponding  to  a new  incident  angle  and  azimuth  at  the  source  and 
tracing  it  through  the  plane  parallel  model  with  a Tocal'  ray  parameter 
determined  by  the  station  azimuth  and  distance  and  structure  dip. 
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An  eastward  dipping  model  is  shown  because  of  the  suggestions  made 

by  McKenzie  and  Julian  (1971)  and  Lin  (1974)  for  an  eastward  dipping  slab. 

Figure  7 shows  the  first  18  seconds  of  eight  representative  P waves 
of  Figure  2.  The  effect  of  layer  dip  manifests  itself  primarily  in  those 
very  stations  which  were  important  in  resolving  the  LVZ,  e.g.,  KEV,  1ST, 

COP,  and  NNA.  Stations  in  which  direct  P is  predominant,  such  as  ATL 
and  GUA,  show  little  change  with  increasing  dip.  The  main  effect  of 
layer  dip  for  all  eastern  azimuths  is  to  have  up-going  rays  start  from 
the  focal  sphere  in  a more  easterly  direction.  Differences  between  the 
focal  area  azimuth  and  station  azimuth  can  be  up  to  30! 

This  change  in  the  relative  amplitude  between  down-  and  up-going  rays 

due  to  the  modified  radiation  pattern  coupled  with  changes  in  the  reflection  coefficie 
for  this  particular  source  and  structure  model  produces  the  effect  of  a 
direct  tradeoff  between  layer  dip  (to  the  east)  and  LVZ  interface  contrasts. 

Although  there  is  no  way  of  knowing  from  these  data  whether  planar  dipping 
structure  exists  at  Puget  Sound  it  is  very  evident  that  it  could  have  a 
very  important  effect  on  the  particular  velocity  contrasts  inferred  from 
the  horizontally  layered  model.  A 50%  reduction  could  easily  be  possible, 
say  for  the  6 = 20°  case,  but  the  general  characteristics  of  the  LVZ 
must  still  be  included  in  the  earth  model. 


Modelling  the  rotated  SH  waves  revealed  little  about  the  source 
or  structure.  Figure  5 shows  the  six  observed  SH  waves  with 
corresponding  synthetic  seismograms  computed  by  ray  techniques  using 
the  PS-9  model.  Basically,  the  observed  SH  waves  are  very  simple 
showing  only  the  direct  S.  Most  of  the  "glichiness"  of  these  waves 
is  due  to  digital  noise  in  the  rotation  process.  Since  upgoing  S 
is  relatively  small  for  these  stations  none  of  the  major  discontinuities 
of  PS-9  are  directly  observable  except  for  the  free  surface.  The 
phase  sS  is  not  readily  apparent  in  the  observations  although  it  is 
theoretically  small.  Several  possibilities  can  be  presented.  First, 
the  S-velocity  structure  of  PS-9  may  be  completely  wrong  so  that  sS 
time  is  significantly  different.  However,  examination  of  the 
SH  observations  yields  no  consistent  arrival  at  any  other  time. 

A second  possibility  is  that  local  receiver  crustal  effects  such  as 
S-coupled  PL  waves  may  contaminate  the  tangential  component.  S-waves 
are  notorious  for  this  and  the  position  of  sS  relative  to  S makes 
this  possibility  very  probable  (Helmberger  and  Engen,  1974).  A 
third  and  very  likely  possibility  is  that  the  earth  model  near  the 
source  is  deficient.  It  could  be  that  anelastic  attenuation  plays 
an  important  role.  Some  short  period  results  will  be  presented  below 
which  support  this  speculation. 

Scaling  the  synthetic  direct  P and  SH  waves  to  the  long  period 

data  yields  a moment  determination  for  the  Puget  Sound  event.  A 

2 6 

value  of  1.3  ± 0.6  x 10^  dyne-cm  is  obtained  from  21  P waves  (see 

Figure  3.2).  The  uncertainty  given  is  the  standard  error.  The  six 

2 6 

SH  waves.  Figure  8,  yield  an  average  value  of  1.7  ± 0.6  x 10  dyne-cm, 
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slightly  higher  than  the  P waves  but  well  within  the  uncertainties  of 

the  geometrical  spreading  correction  and  assumptions  on  attenuation.  The 

2 6 

average  moment  including  both  P and  S waves,  is  1.4  ± 0.6  x 10^  dyne-cm. 

The  scatter  in  amplitudes,  yet  not  in  waveshape,  can  be  considerable  even 
for  nearby  stations.  For  example,  compare  1ST  and  ATU,  Figure  2. 

An  attempt  was  made  to  utilize  the  many  short  period  observations 
from  this  event  to  determine  timing  and  amplitude  information  to  help 
pin  down  the  inferred  structural  discontinuities.  Although  the  attempt 
eventually  fell  short  of  its  original  goals  some  interesting  effects  between 
the  direct  and  reflected  phases  were  observed. 
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An  easy  way  to  change  a short  period  record  into  a form 
suitable  for  stacking  is  to  compute  an  envelope  of  the  signal.  Once 
the  envelope  is  found  for  several  records  a suitable  amplitude 
normalization  is  performed  and  the  traces  averaged  with  all  seismo- 
grams lined  up  with  respect  to  their  first  arrivals.  A convenient  way 
to  compute  the  type  of  envelope  needed  here  is  to  take  the 
instantaneous  amplitude  of  the  analytic  signal  (Fambach,  1975). 

The  analytic  signal  is  defined  by  (Bracewell , 1965) 


S(t)  = f(t)  - iFH.(t) 


= |s(t)| 


eia(t) 


(1) 


where , 

f(t)  = observed  time  series 

S(t)  = analytic  signal 

F„.(t)  = Hilbert  transform  of  f(t) 

Hi 

a(t)  = time  varying  phase 
| S ( t ) | = instantaneous  amplitude. 

The  Hilbert  Transform  in  equation  (1)  can  simply  be  thought  of  as 
a convolution  of  -1/irt  with  f(t).  Figure  9 shows  an  example  of 
taking  the  instantaneous  amplitude,  |S(t)|,  of  a short  period  vertical 
P wave  recorded  at  AFI . This  method  preserves  the  times  and  relative 
amplitudes  of  the  arrivals  and  even  makes  the  record  a little  easier 
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to  visually  interpret.  By  no  means  does  it  purport  to  add  any 
information;  it  simply  makes  the  sesimograms  easier  to  work  with. 

The  normalization  was  done  by  setting  the  area  of  the  first  10  seconds 
of  the  envelope  to  unity.  This  tended  to  boost  the  relative  amplitude 
of  sharp  arrivals  to  that  of  emergent  arrivals. 

Figure  10  - shows  the  results  obtained  from  a stacking  experiment. 

The  traces  labelled  'P'  are  stacked  envelopes  of  short  period  vertical 

P waves  grouped  as  a function  of  range.  The  range  interval  is 

indicated  under  each  group  name  and  the  number  in  parentheses 

represents  the  number  of  seismograms  used  per  group.  Table  1 

indicates  the  particular  stations  used  in  this  grouping.  The  P wave 

stacked  traces  show  several  interesting  effects.  In  the  first  10  sec 

two  prominent  arrivals  are  usually  evident  and  are  designated  P and 

p^P  in  the  diagram.  The  first  arrival  is  interpreted  to  be  the 

direct  P wave  with  p P being  the  reflection  from  the  bottom  of  the 

m 

LVZ.  The  timing  agrees  with  the  long  period  result.  This  is  not 
too  useful  except  to  show  that  the  reflection  has  approximately  the 
same  frequency  content  as  the  direct  wave.  This  implies  that  the 
contrast  at  this  boundary  must  be  sharp,  probably  less  than  2 km  in 
transition.  The  most  striking  effect  these  traces  demonstrate,  however, 
is  the  very  conspicuous  absence  of  pP , especially  at  the  smaller 
epicentral  distances.  A simple  glance  at  the  long  period  waveforms 
in  Figure  2 reveals  that,  except  for  the  western  stations,  pP 
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is  at  least  as  big  as  direct  P and  usually  two  or  three  times 
bigger.  The  relative  amplitude  of  pP  in  groups  a,  b,  and  c of 
Figure  10  is  at  least  half  that  of  direct  P.  Not  until  stations 
with  ranges  greater  than  75°  are  considered  does  pP  start  becoming 
apparent.  This  effect  can  be  inferred  to  come  from  earth  structure 
by  using  the  information  gleaned  from  p^P.  Because  the.  frequency 
content  of  upgoing  rays  is  similar  to  that  of  the  downgoing  rays 
a source  effect,  such  as  directivity,  can  automatically  be  eliminated. 
All  reflections  should  be  apparent  unless  structure  or  attenuation 
effects  cancel  the  arrivals. 

An  explanation  for  this  observation  can  be  found  in  the  geologic 
and  structural  framework  of  Puget  Sound.  Figure  11  is  a map  of 
southern  Puget  Sound  with  the  epicenter  plotted  near  the  center 
(after  Algermissen  and  Harding,  1965).  The  concentric  circles  arc 
contours  of  the  position  pP  hits  the  surface  for  a particular 
station  distance.  These  were  computed  from  ray  theory  for  model 
PS-9  using  the  expression 

xr  = U h~  Thi  (2> 

j=i  M J 

where 

x = radius  from  epicenter 

r 1 

p = ray  parameter  considered 
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Th.  = thickness  of  the  jtl  layer 

The  summation  is  for  all  layers  between  the  hypocenter  and  the  free 
surface.  The  hachured  lines  across  Figure  11  represent  a simplified 
version  of  structural  discontinuities  from  the  gravity  and  magnetic 
interpretation  of  Danes  ej^  aJL  (1965).  They  are  dashed  where 
questionable.  The  tick  marks  are  on  the  down-thrown  side  of  the 
inferred  fault.  The  large  numbers  in  parentheses  are  inferred  depths 
to  crystalline  basement  from  the  surface.  This  map  readily 
demonstrates  the  great  structural  relief  and  complexity  in  the 
near-surface  geology.  The  contours  of  pP  are  only  approximate  when 
considered  from  this  approach  since  they  were  calculated  using  a 
plane  layered  model,  an  assumption  which  clearly  breaks  down  for  the 
near  surface  layers.  Nevertheless,  they  should  be  useful  to  first 
order.  Comparing  the  distance  ranges  and  the  structure  gives  a 
possible  explanation  for  the  amplitude  behavior  of  short  period  pP. 

For  all  stations  at  distances  less  than  75°  pP  bounces  at  points 
which  have  thick  sections  of  sediment.  For  stations  with  ranges 
greater  than  that,  pP  bounces  within  the  boundaries  of  the  large 
horst  or  near  its  edge  where  sediment  thicknesses  are  presumably 
smaller.  This  observation  seems  to  correlate  with  the  effects  seen 
in  the  traces  of  Figure  10.  Possibly,  then,  the  attenuation  of  pP 


relative  to  P comes  from  either  structurally  disturbed  sediments 
causing  wave  scattering  or  perhaps  from  anelastic  attenuation. 
Either  mechanism  is  feasible.  An  effective.  ’ Q * can  be  computed 
from  these  observations  by  using  equation  (3).  A conservative 

estimate  of  the  amplitude  attenuation  factor,  a,  is  about  1/3  for 
a .75  sec  short  period  sinusoid,  and  with  a total  travel  time,  T, 
of  17  sec  gives 


Q 


eff 


-ti)T 

2 In  a 


(3) 


= 65 


This  amounts  to  convolving  pP  with  another  Q operator  on  top  of  the 
first  (T/Q  = 1.0)  with  a travel-time  to  Q ratio  of  0.3.  This 
particular  Q operator  would  have  negligible  effect  on  a long  period  P 
waveform. 

The  realization  that  Puget  Sound  is  structurally  heterogeneous 
can  explain  why  the  shape  details  of  pP  are  not  fit  by  PS-9  for  all 
azimuths  (see  Figure  2).  In  fact,  it  is  surprising  that  the  model 
works  at  all.  Perhaps  the  pathological  case  of  GPU , Figure  2, 
can  have  an  explanation  in  this  light. 

Returning  to  Figure  10  one  last  point  can  be  made  using  the 
short  period  stacking  procedure.  The  bottom  trace,  labeled  'S', 
is  the  average  of  both  horizontal  short  period  S wave  components 


at  seven  stations  (Tab! 


1).  The  result  of  the  averaging  revealed 
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two  arrivals  roughly  3 sec  apart.  The  first  is  interpreted  to  be 
direct  S.  Comparisons  with  each  corresponding  long  period  component 
revealed  that  the  short  period  S started  after  any  long  period  S-P 
precursors  and  correlated  with  where  the  long  period  arrival 
was  sharpest.  The  second  arrival,  because 

of  its  timing,  is  interpreted  to  be  the  reflection  from  the  bottom 
of  the  LVZ.  This  arrival  did  not  occur  at  all  stations  and  there 
were  significant  waveform  variations  between  stations  to  make 
this  a tentative  conclusion.  It  does  give  weak  evidence  against 
multiple  source  complications  through  the  agreement  in  travel  times 
for  the  inferred  reflected  S phase. 

Discussion  of  the  earthquake  results 

The  source  model  presented  here  is  an  exceedingly  simple  one. 

The  orientation  agrees  with  previous  studies  and  the  time  function 
is  short  and  uncomplicated.  It  is  certainly  quite  possible  that 
there  are  other  source  complications  but  an  earth  model  has  been 
presented  which,  if  anything,  demonstrates  the  problems  one  has  to 
deal  with  in  order  to  discern  these  effects. 

The  earth  model,  on  the  other  hand,  is  relatively  -omplicated. 

In  some  respects  it  is  heartening  to  see  how  some  parts  of  the  model 
agree  with  the  geology  and  previous  geophysics.  For  example,  the 
inferred  crustal  thickness  is  about  15  Kjn.  The  refraction 
results  are  consistent  with  this  including  a recent  study  using 
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local  mine  blasts  recorded  by  the  Puget  Sound  seismic  array.  This 
particular  study  also  obtained  a shallow  depth  to  Moho  of  15  to  20  km 
(H.  Zuercher,  1975).  The  low  velocity  layers 

at  the  top  of  the  model,  needed  to  widen  pP  and  attenuate  sP,  correlate 
nicely  with  the  thick  sections  of  presumably  Tertiary  sediments  under  the 
blanket  of  glacial  till. 

On  the  other  hand,  discovery  of  the  massive  LVZ  in  the  uppermost  mantle 
by  these  means  presents  many  problems  in  uniqueness,  although  the  model 
produces  a good  fit  to  all  the  P waves.  Clearly,  a check  on  this  conclusion 
is  desirable.  One  way  to  do  this  would  be  to  examine  another  earthquake 
to  see  if  the  same  structure  effects  are  observable.  The  1949  event  would 
be  the  logical  choice  since  it  occurred  nearby  and  at  a depth  of  70  km. 

The  records,  however,  are  not  easily  available  and  the  orientation  of  the 
event,  relative  to  the  station  coverage,  might  not  he  appropriate  for 
observing  these  effects.  Because  the  inferred  structural  boundaries  of 
the  LVZ  are  so  distinct  they  should  be  directly  observable  from  upcoming 
phase  conversions  and  reverberations  using  teleseismic  sources  for  stations 
situated  over  the  structure.  Fortunately,  there  are  several  long  period 
WWSS  and  Canadian  network  stations  nearby  to  test  the  model.  I'he  WVJSSN 
station  at  Corvallis,  Oregon,  was  studied  in  this  manner  (Langston,  1976, 
this  volume)  and  although  the  details  of  the  earth  mode!  are  diffen  al  fro-: 
PS-9,  the  distinct  LVZ  was  again  found. 

It  is  interesting  to  note  that  in  a recent  set  of  papers  on  earthquake 
travel  times  in  Puget  Sound,  Crosson  (1976)  found  a possible  LVZ  although 
he  states  it  may  not  be  resolvable  from  his  data.  Further  discussion 
will  be  given  in  Langston  (1976  b,  this  volume). 
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Conclusions 

A major,  distinct  low  velocity  2one  is  inferred  in  the  uppermost 
mantle  under  Puget  Sound  by  modelling  reflected  arrivals  in  the  long 
period  P waveforms  from  the  1965  Puget  Sound  earthquake.  The  LVZ  occurs 
in  the  depth  interval  of  41  to  56  km  and  is  sandwiched  between  layers 
having  mantle  velocities.  The  velocity  contrast  at  the  bottom  interface 
of  the  LVZ  must  be  on  the  order  of  2.5  km/sec  as  inferred  from  the 
theoretical  relative  P amplitudes  from  the  point  dislocation  source, 
assuming  horizontal  structure.  If  eastward  dipping  structure  is  allowed, 
the  velocity  contrasts  in  the  LVZ  trade-off  directly  with  dip. 

The  crust  of  southern  Puget  Sound  is  thin,  less  than  15  km,  and 
has  a thick  sediment  section  near  the  surface  as  inferred  from  the  chapes 
of  the  pP  and  sP  phases.  This  is  consistent  with  other  geophysical 
studies. 

The  earthquake  source  parameters  determined  by  the  P and  SH  waveform 
study  are  the  following:  70°  dip  to  the  east;  344°  strike;  -75°  rake, 

making  it  a normal  dip-slip  fault,  consistent  with  previous  studies; 

2 () 

seismic  moment  of  1.4  ± 0.6  x 10  dyne-cm;  63  km  depth;  and  a simple 
time  history  represented  by  a triangular  far-field  source  function  with 
a rise  time  of  0.5  sec  and  a fall  off  of  2.5  sec.  The  depth  is  constrained 
by  pP  and  sP  travel  times. 

The  synthetic  P and  SH  waveforms  computed  from  the  given  sour'  .v.i.l 
earth  models  reproduce  the  major  azimuthal  variations  which  occur  in  the 
data  waveforms. 

A short  period  stacking  procedure  utilizing  the  instantaneous 
amplitude  of  the  analytic  signal  is  presented  and  used  to  find  reflected 
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arrivals.  An  apparent  attenuation  is  found  for  short  period  versus  long 
period  pP  and  is  correlated  with  sediment  thickness  at  the  reflection 
point.  The  amplitude  discrepancy  yields  an  effective  Q of  about  65 
for  the  earth  structure  above  the  earthquake. 
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Table  1 

Station  Information 


Station 

A(0) 

AZ(0) 

BAZ(0) 

Components 

AFI 

A 

75.3 

229.7 

32.3 

P,  P(d) 

AKU 

53.2 

30.1 

304.7 

P(b) 

ANP 

87.8 

305.5 

37.6 

P 

ANT 

84.8 

133.6 

327.5 

P,  p(d) 

ARE 

78.2 

130.5 

327.5 

s 

ATL 

31.8 

102.2 

307.4 

p,  p(a) 

ATU 

89.8 

26.2 

337.7 

P,  p(d) 

BEC 

45.7 

87.6 

306.7 

P,  p(b) 

CAR 

59.2 

109.5 

319.4 

P,  p(b),  s 

COP 

70.5 

25.3 

329.2 

P,  S 

ESK 

65.6 

33.3 

319.3 

P (c) 

CDH 

39.4 

31.6 

272.9 

P,  p (a) , S 

CEO 

33.  7 

87.8 

299.5 

P,  P (a) 

GTE 

55.5 

140.0 

334.1 

P (b) 

CUA 

82.0 

281.2 

43.3 

P 

IINR 

88.7 

254.  7 

41.6 

P,  p(d) 

1ST 

88.3 

21.  3 

340.9 

P,  P(d) 

KEV 

61.0 

11.7 

336.6 

p,  r(e) 

KON 

66.3 

24.4 

326.5 

p 

KTC 

49 . 5 

25.9 

298.4 

P,  p (b) , s 

J.PB 


SO . 0 


127. 8 


026.0 


I’ , p (cl)  , n 


Tab  Lc  1 


(continued) 


Station 

A(0) 

AZ(0) 

BAZ^0) 

Components 

MAL 

79.5 

46.2 

322.4 

P 

NNA 

71.7 

132.7 

329.3 

P,  s 

NOR 

45.7 

11.4 

293.8 

s 

NUR 

69.2 

16.  8 

336.6 

P,  P (c)  , S 

OGD 

34.3 

82.8 

297.0 

p (a) 

PDA 

67.4 

58.6 

313.0 

p 

PTO 

74.1 

46.0 

319.7 

p,  p(c),  s 

QUI 

60.8 

127.5 

327.4 

P 

SCP 

32.2 

85.2 

296.9 

P 

SHA 

31.0 

110.2 

312.3 

P 

SJG 

54. 1 

102.8 

315.8 

P 

STU 

75.4 

30.  8 

328.3 

P,  p(d) 

TOL 

77.2 

44.0 

322.1 

P,  s 

TRN 

62 . 6 

104.  7 

318.  1 

p.  p(c) 

UME 

65. 3 

17.  3 

332.9 

p(c),  S 

VAL 

65.3 

39.2 

316.1 

p(c)  , S 

WES 

35.9 

78.8 

295.9 

P,  P (a) 

long  period 

vertical  P 

waveform  or 

first  mot 

ion 

: short  period  vertical 

P waveform 

Letter 

in  parenthesis 

designates  stacking 

group . 

long  period 

tangent i al 

S waveform 

short  period 

hori  ;:onta  1 

S waveform' 

used  in 

si aek i ng 

Table  1 


Structure  Models 


3)  Layer  Over  Halfspace 


Layer  # 

V 

P 

V 

s 

P 

Th 

1 

3.0 

1.2 

1.8 

2.0 

2 

4.0 

2.0 

2.0 

2.0 

3 

5.0 

2.5 

2.2 

1.0 

4 

6.0 

3.0 

2.  3 

1.0 

5 

7.0 

3.5 

2.4 

1.0 

6 

00 

r*. 

4.0 

3.2 

28.0 

7 

6.5 

3.0 

2.7 

10.0 

8 

5.5 

2.7 

2.6 

10.0 

9 

8.0 

4.3 

3.2 

- 

1 

2.0 

1.0 

2.5 

1.2 

2 

3.0 

1.5 

2.5 

1.2 

3 

4.5 

2.3 

2.5 

2.5 

4 

6.8 

3.9 

2.9 

4.0 

5 

7.4 

4.2 

3.0 

4.5 

6 

7.8 

4.3 

3.2 

27.  3 

7 

6.5 

3. 1 

2.9 

7.0 

8 

5.5 

2.9 

2.7 

8.0 

9 

8.0 

4.6 

3.  2 

- 

1 

6.0 

3.5 

2.7 

10.0 

2 

8.0 

4.6 

3.2 
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Figure  Captions 


Figure  1. 


Figure  2. 


Figure  3. 
Figure  4. 

Figure  5. 
Figure  6. 
Figure  7. 


Index  map  of  western  Washington  and  Oregon  showing  the  epicenter 
of  the  1965  Puget  Sound  event  and  WWSSN  stations  used  in 
the  receiver  structure  determination.  The  dashed  line  is 
the  approximate  extent  of  the  Tertiary  eugeosync I ine , after 
Snavely  and  Wagner  (1963). 

Comparison  of  synthetic  and  observed  P waveforms  for  the  final 

source  and  earth  structure  models.  The  observed  is  displayed 

directly  above  the  synthetic  at  each  station  with  the  calculated 

26 

moment  indicated  next  to  the  traces  (x  10  dyne-cm) . The  focal 
plot  is  for  the  bottom  hemisphere. 

PS-9  earth  model  for  Puget  Sound. 

Comparison  of  P wave  synthetic  seismograms  at  five  representative 
stations  for  the  one  layer  model  (top)  and  PS-9  (bottom).  The 
phase  identifications  are  referred  to  in  the  text. 

KEV  ray  summation  showing  the  effects  caused  by  ref  1 e'-r  ions  from 
interfaces  hounding  the  inferred  LVZ. 

Comparison  of  P wave  synthetic  seismograms  computed  using  tin  PS-1 
and  PS-9  earth  models  (Table  2)  with  the  ohserv.it  ion  it  KFV  (middle). 
Parameter  study  in  which  the  horizontal  layer  interfaces  at  the  IV 
of  PS-9  are  allowed  to  dip.  At  each  slat  ion  if*  to'  mm.,  tin 
observed  P wave  and  just  below  it  is  tin  svntln  ’ foi  P:  9. 

Below  this  waveform  are  two  synthetics  tor  i 1 1 » 1-  and  1 I 


dipping  LVZ. 
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Figure  8. 


Figure  9. 


Figure  10. 


Figure  II. 


Comparison  of  the  observed  and  final  synthetic  SH  waves.  Same 
scheme  as  in  Figure  2. 

Example  of  finding  the  envelope  of  a short  period  record  using  the 
instantaneous  amplitude  of  the  analytic  signal. 

Stacked  envelopes  of  short  period  P and  S waves.  The  arrows  point 
directly  to  the  phase  being  identified  except  for  the  pP  arrow, 
which  designates  the  theoretical  arrival  time. 

Map  of  the  southern  Puget  Sound  area  showing  the  epicenter  for  the 
1965  event  (after  Algermissen  and  Harding,  1965).  The  concentric 
circles  are  contours  of  pP  reflection  points  as  a function  of  epicentral 
distance.  The  hachured  lines  are  faults  inferred  by  gravity  and 
magnetics  from  Danes  et  al . (1965)  with  the  tick  marks  on  the 
downthrown  side.  The  numbers  in  parenthesis  are  depth  t basement  in  km. 
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ABSTRACT 

A dataset  of  forty-one  moderate  and  large  earthquakes  has  been 
used  to  derive  scaling  rules  for  kinematic  fault  parameters.  If  effec- 
tive stress  and  static  stress  drop  are  equal,  then  fault  rise  time,  x, 

1/2  3/2 

and  fault  area,  S,  are  related  by  t = 16S  /(7n  (3),  where  0 is  shear 

velocity.  Fault  length  (parallel  to  strike)  and  width  (parallel  to  dip) 

are  empirically  related  by  L = 2W.  Scatter  for  both  scaling  rules  is 

about  a factor  of  two.  These  scaling  laws  combine  to  give  width  and 

rise  time  in  terms  of  fault  length.  Length  is  then  used  as  the  sole 

free  parameter  in  a Haskell  type  fault  model  to  derive  scaling  laws 

relating  seismic  moment  to  Mg  (20  sec.  surface  wave  magnitude),  Mg 

to  S and  (1  sec  body  wave  magnitude)  to  M^.  Observed  data  agree 

well  with  the  predicted  scaling  relation.  The  "source  spectrum"  depends 

on  both  azimuth  and  apparent  velocity  of  the  phase  or  mode,  so  there  is 

a different  "source  spectrum"  for  each  mode,  rather  than  a single  spectrum 

for  all  modes.  Furthermore,  fault  width  (i.e.  the  two  dimensionality 

of  faults)  must  not  be  neglected.  Inclusion  of  width  leads  to  different 

average  source  spectra  for  surface  waves  and  body  waves.  These  spectra 

in  turn  imply  that  and  Mg  reach  maximum  values  regardless  of  further 

increases  in  L and  seismic  moment.  The  relation  from  this  study 

differs  significantly  from  the  Gutenberg-Richter  relation,  because  the 

G-R  equation  was  derived  for  body  waves  with  a predominant  period  of  about 

3 sec  and  thus  does  not  apply  to  modern  1 sec  m determinations.  Previous 
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investigators  who  assumed  that  the  G-R  relation  was  derived  from  1 sec 

data  were  in  error.  Finally  averaging  reported  rupture  velocities 

yields  the  relation  V = .72R. 

K 

INTRODUCTION 

The  purpose  of  this  paper  is  to  examine  empirical  relations  between 
gross  fault  parameters  and  the  agreement  of  these  relations  with  theo- 
retical models  of  seismic  sources.  The  gross  parameters  to  be  studied 
are  fault  length,  width  and  rise  time,  rupture  velocity  m^  and  Mg,  and 
seismic  moment.  Data  from  other  investigators’  studies  of  individual 
earthquakes  are  used  to  study  scaling  of  source  parameters  in  an  approxi- 
mate way.  In  general  the  data  are  consistent  with  fault  width  scaling 
proportionately  to  fault  length  and  rise  time  scaling  proportionately 
to  the  square  root  of  fault  area.  This  scaling  can  then  be  used  to  find 

m.  - M , log  M - M and  log  S - ” relations.  Some,  of  those  relations 
o s o s s 

have  been  studied  by  Kanamori  and  Anderson  (1975b), who  provided  a theo- 
retical basis  for  many  of  the  empirical  relationships  used  in  seismology. 

Tsuboi  (1956)  was  the  first  investigator  to  utilize  similarity,  the 
concept  of  relating  earthquakes  of  different  sizes  by  a one  parameter 
model.  By  assuming  that  the  horizontal  dimensions  of  the  earthquake 

source  volume  were  three  times  larger  than  the  vertical  dimension, 

2 15 

Tsuboi  derived  from  the  relation  E = pc  A / 6 , where  F.  is  released  energy, 
M is  average  rigidity,  c is  average  strain  drop  and  A is  aftershock 
area.  Such  approximate  scaling  relations,  as  first  pointed  out  bv 
Tsuboi,  require  that  the  physics  of  material  failure  be  identical  for 


large  and  small  earthquakes.  If  that  assumption  is  generally  true  and 
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if  earthquakes  tend  to  be  georietrically  similar,  then  it  follows  that 
fault  length  and  width,  and  final  dislocation  all  will  scale  together. 
Differences  in  material  properties  will  weaken  the  exactness  of  the 
similarity  when  earthquakes  from  two  different  regions  are  compared, 
but,  in  an  approximate  sense,  similarity,  as  is  shown  by  the  data  pre- 
sented below,  is  a valid  concept. 

The  first  paper  to  systematically  relate  observed  gross  seismic 
source  parameters  to  the  source  spectrum  was  the  now  classic  work  of 
Aki  (1967).  Although  the  results  presented  in  this  paper  modify  his 
results,  the  methodology  and  basic  outlook  are  similar  to  Aki’s.  Later 
Brune  (1970,  1971)  contributed  to  the  understanding  of  seismic  source 
spectra. 

Similarity  between  earthquakes  is  a dynamic  as  well  as  a static 
concept.  Mot  only  the  final  static  parameters,  but  also  the  spectral 
shape  of  the  equivalent  source  time  function,  scales  with  fault  length. 
Spectral  similarity  can  best  be  demonstrated  by  comparing  two  earth- 
quakes with  identical  location  and  focal  mechanism,  but  different 
magnitude.  Such  a comparison  ensures  that  seismograms  from  both  events 
will  be  affected  equally  by  the  medium  response,  so  that  all  differ- 
ences between  the  records  will  be  from  the  source  effects. 

Observational  studies  of  similar  pairs  of  earthquakes  have  been 

made  by  Berckhemer  (1962)  whose  results  were  interpreted  by  Aki  (1967) 
_2 

to  support  Aki's  w model.  Tsujiura  (1973)  studied  groups  of  events 
from  various  regions,  concluding  that  most  data  were  In  accord  with 
Aki's  w model,  but  that  some  were  better  fit  by  an  u>  ^ or  tiP  model. 


One  cannot  directly  compare  spectral  characteristics  of  source 
mechanisms  from  different  regions  without  first  correcting  the  seismo- 
grams for  transmission  effects.  Removing  the  effects  of  medium  response 
will  usually  require  use  of  synthetic  seismogram  methods.  We  assume 
however  that  one  can  compare  logarithmic  fault  parameters,  such  as  , 

Mg  or  log  L,  for  events  in  different  regions.  These  comparison  are 
made  with  the  intention  of  looking  at  order  of  magnitude  relationships 
rather  than  details. 

In  this  paper  we  will  look,  at  scaling  relations  between  five  sets 
of  such  logarithmic  parameters:  log  L vs.  log  W (fault  length  vs. 

width),  log  t vs.  log  S (fault  rise  time  vs.  area),  Ms  vs.  log  M0 
(surface  wave  magnitude  vs.  seismic  moment),  1%  vs.  Ms  and  log  S vs.  Ms. 

What  will  be  shown  are  not  exact  correlations,  but  rather  trends  which 
appear  to  be  applicable  to  most  earthquakes.  Agreement  between  the 
simple  model  used  in  this  paper  and  the  data  are  quite  good.  (Long 
narrow  transform  faults,  such  as  the  San  Andreas,  are  1 separate  class  of 
faults  which  are  not  considered  in  this  paper.) 

OBSERVATIONAL  DATASET 

The  earthquake  data  shown  in  Table  1 are  from  the  s.>:\  forty-one  shallow 
events  used  by  Kanamori  and  Anderson  (1975b).  All  values  lor  M are 
from  their  paper;  the  sources  for  all  other  observational  parameters 
are  given  in  the  table  of  references.  Except  for  minor  differences 
which  are  primarily  due  to  the  use  of  slightly  different  references, 
data  for  M_  and  S are  equivalent  to  Kanamori  and  Andersens '.  Each 
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numbered  entry  in  the  table  of  references  corresponds  to  the  earthquake 
with  the  same  number  in  Table  1.  All  but  two  of  the  columns  are 
observational  data;t*,  predicted  rise  time,  and  A a,  calculated  stress 
drop,  will  be  discussed  below.  Length  and  width  have  been  taken  from 
the  references,  or  in  some  cases  estimated.  Length  always  refers  to 
length  along  the  strike,  regardless  of  focal  mechanism;  width  refers 
to  width  along  the  dip.  Average  dislocation  comes  either  from  field 
measurements  or  from  divi  ling  the  moment  (determined  from  seismograms) 
by  the  area  and  on  assumed  value  of  the  shear  modulus. 

For  all  events  since  August  1963  the  value  is  either  taken 
directly  from  the  PDF.  Monthly  Summary , or  calculated  from  the  data  in 
Earthquake  Data  Reports.  As  reported  by  L.  M.  Murphy  in  Bath  (1969)  , 
USCCS  (later  NOAA  and  now  USGS)  asks  for  the  amplitude  of  the  largest 
pulse  (with  period  less  than  three  seconds)  in  the  first  five  cycles  of 
the  teleseismic  P or  Pn  arrival.  The  values  of  A and  T are  then  used 
in  the  Gutenberg-Richter  formula 

= log10(A/T)  + Q,  (1) 

to  derive  for  each  station.  Values  more  than  0.7  magnitude  units 
from  the  mean  are  deleted  and  the  final  average  is  then  taken. 

Estimates  of  rise  times  typically  were  made  by  fitting  the  first 


upswing  on  long-period  local  records  to  synthet ic  seismograms  calcu- 
lated using  the  Haskell  (1969)  whole  space  model  at  one  or  two  stations. 
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Clearly  it  would  be  desirable  to  use  synthetics  made  for  more  realistic 
models  of  earth  structure,  but  they  have  not  yet  been  calculated  for 
these  events.  Uncertainties  due  to  the  tradeoff  between  rise  time  and 
rupture  velocity  and  due  to  the  model  may  combine  to  cause  errors  which 
cannot  be  estimated.  In  some  cases,  such  as  the  Tot  tori  earthquake 
(Kanamori,  1972b),  rupture  velocity  and  rise  time  are  independently 
constrained. 


LENGTH  VERSUS  WIDTH 

Fault  length  is  plotted  against  fault  width  in  Figure  1.  It  can 
be  seen  that  (with  considerable  scatter)  observational  data  demon- 
strate that  L = 2W.  In  Figure  1,  the  numbers  refer  to  earthquakes  in 
Table  1.  Intraplate  events  are  plotted  as  open  circles  and  interplate 
events  as  solid  circles.  (This  convention  is  also  used  in  all  later 
figures.)  There  is  not  any  clear  difference  between  the  interplace  and 
intraplate  groups.  Abe  (1976)  has  independently  found  I.  = 2W  for  a 
dataset  of  Japanese  earthquakes. 

The  Haskell  model  uses  L as  the  direction  in  which  rupture 
propagates,  while  L was  measured  along  the  direction  of  the  strike  for 
Figure  1.  It  is  implicitly  assumed,  then,  that  for  these  41  events 
rupture  propagated  parallel  to  the  strike.  This  almost  certainly  false 
for  some  thrust  events  such  as  San  Fernando  ( doom  and  F.oback  (1974)  » 
Trifunac  (1974),  J'ikumo  (1971b) ) ,and  may  well  be  lalso  for  events  like 
Hemuro-Oki  L < W.  in  spite  of  these  exceptions  it  seems  that  rupture 
usually  propagates  parallel  to  the  strike,  especially  for  strike-slip 


fan! ts. 
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RISE  TIME  VERSUS  THEORETICAL  PREDICTIONS 
Kanamori  (1972b)  showed  that 

• _ 5 n 6 (2) 

D = — E;  0 — 

T e,o  p 

where  D is  dislocation  velocity,  5 is  average  dislocation,  t is  rise 

time,  and  o is  effective  stress.  If  one  assumes  that  effective 
e,o 

dynamic  stress  is  equal  to  static  stress  drop,  Ao,  this  assumption  can 
be  tested  by  comparing  observed  rise  times  to  the  theoretically  pre- 
dicted rise  time 


T* 


pD 

BAo 


(3) 


One  can  obtain  stress  drop  in  closed  form  for  only  a few  simple 
models.  The  most  straightforward  of  these  is  the  circular  crack  with 
constant  stress  drop  discussed  by  Keiles-Borok  (1959).  For  that  model 
stress  drop  is  given  by 

Ao  = 7tt3/2pD/(16»'’s)  = 7Mo/(16(LW/it)3/2)  (4) 


Although  this  formula  does  not  give  the  exact  stress  drop  for  the 
rectangular  fault  model,  it  follows  from  the  work  of  Sato  (1972)  that 
this  is  a good  approximation.  If  we  substitute  (4)  in  (3),  where  S is 
the  area  of  the  rectangular  fault  and  D is  average  dislocation,  we 


obtain 
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t*  = 16S1/2/(7h3/23)  (5) 

The  values  of  A a and  t*  in  Table  1 were  calculated  using  (4)  and  (5) , 
respectively. 

Figure  2 is  a plot  of  observed  versus  predicted  (from  (5))  rise 
times  for  a number  of  earthquakes.  It  can  be  seen  that,  again  with 
considerable  scatter,  observational  and  theoretical  rise  times  are  in 
agreement.  Abe  (1975b)  reached  a similar  conclusion  from  a dataset  of 
five  Japanese  earthquakes. 

The  agreement  between  theoretical  and  observed  rise  times  has 
important  implications  for  engineering  seismology.  The  only  observa- 
tional parameter  required  in  (5)  is  fault  area,  which  frequently  can  be 
estimated  from  geological  data.  If  total  dislocation  can  also  be 
determined  from  geological  field  work,  then  particle  velocity  near  the 
fault,  an  important  parameter  in  engineering  seismology,  can  be  reliably 
estimated.  This  is  potentially  of  great  values  in  areas  lacking  in 
historical  seismicity  or  good  instrumental  data. 


RUPTURE  VELOCITY 

Table  1 lists  rupture  velocities  reported  by  various  investigators. 
These  values  were  determined  from  matching  synthetic  seismograms  to 
local  records  or  from  surface  wave  analysis.  To  a certain  extent  then, 
these  values  are  model  dependent.  Some  also  may  be  affected  by  the 
difficulty  in  resolution  between  rise  time  and  rupture  velocity.  Never- 
theless, these  measurements  probabLy  represent  a good  average  sample  of 
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rupture  velocity  measurements.  If  one  picks  values  of  3 ranging  from 
3.5  km/sec  for  shallow  crustal  events  to  A. 5 km/sec  for  events  breaking 
the  entire  lithosphere,  one  then  can  calculate  that  the  average  value  of 
(Vr/3)  is  0.72. 

Archuleta  and  Brune  (1975)  found  V /3  = 0.7  in  experiments  on 

R 

fracture  of  prestressed  foam  rubber.  Their  measured  value  was  for  the 
surface  of  the  foam  rubber,  but  if  one  assumes  infinite  rupture  velocity 
along  the  dip,  their  result  agrees  very  well  with  the  result  VR/3  = 0.72 
observed  for  earthquakes.  (Their  minimum  possible  value  for  VR/3  at 
depth  is  0.636.)  Agreement  between  the  earthquake  and  foam  rubber 
rupture  velocities  may  be  fortuitous  or  may  be  caused  by  a common 
physical  friction  mechanism. 

CHOICE  OF  FAULT  MODELS 

All  "deterministic"  source  models  specify  some  (nearly  always 
kinematic)  conditions  at  the  source,  which  then  fix  via  the  repre- 
sentation theorem  of  de  Hoop  (1958)  and  Burridge  and  Knopoff  (1964),  the 
complete  time  history  at  every  point  in  the  medium.  ( Aki  (1967,  1972) 
and  Haskell  (1966)  proposed  "statistical"  models  in  which  only  the 
amplitude  spectrum  at  the  source  function  is  specified.  Since  we  will 
be  looking  at  dislocation  rise  times,  these  statistical  models  are  not 
appropriate  choices.)  Typically  the  source  theory  papers  calculate 
seismograms  for  an  isotropic  homogeneous  whole  space.  Since  our  inter- 
est is  in  logarithmic  source  parameters,  we  will  assume  the  whole  space 


models  are  adequate. 
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In  most  deterministic  source  models,  either  fault  dislocation 
(e.g.,  Haskell  (1969),  Mikumo  (1973b))  or  stress  drop  (and  therefore 
fault  dislocation,  e.g.,  Burridge  and  Willis  (1969),  Richards  (1973), 

Sato  and  Hirasawa  (1973)), is  specified,  which  in  turn  gives  displacement 
at  other  points  in  the  medium.  Other  authors,  e.g. , Hanson  et^  al. 

(1974)  and  Andrews  (1975), have  studied  numerical  models  with  friction 
between  the  fault  surfaces. 

All  of  these  models  predict  far-field  pulses  which  scale  linearly 
with  fault  dimensions.  Also  they  all  yield  flat  spectra  at  low  fre- 
quencies and  w n high  frequency  asymptotes  (n  5 2).  Thus  all  of  the 
models  have  at  least  one  "comer  frequency"  (and  some  have  several). 

For  these  models  the  static  or  low  frequency  level,  which  is  proportional 

3 

to  seismic  moment,  grows  as  L . 

We  will  continue  to  use  the  Haskell  (1964,  1969)  model  of  a 
rectangular  fault  (shown  in  Fig.  3)  in  this  paper.  Most  studies  have 
used  this  model  in  the  determination  of  rise  times  from  local  seismo- 
grams. The  basic  Haskell  model  is  a fault  with  length  L,  width  W,  rise 
time  t (linear  ramp  time  function),  final  dislocation  D and  rupture 

velocity  V . Rupture  is  instantaneous  in  the  width  direction  and 
K 

propagates  (starting  at  one  end)  along  the  length  with  velocity  V . 

Some  investigators  have  made  the  natural  extension  to  bilateral  rupture 
propagation. 

Haskell's  (1964)  expressions  are  for  a "one-dimensional"  model  in 
which  width  is  included  only  as  a weighting  factor  in  the  moment. 

Hirasawa  and  Stauder  (1965)  and  Mikumo  (1969)  included  the  complete 
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effect  of  the  width  to  obtain  an  expression  for  spectral  source 
amplitude. 


U (u) 

c 


m„r 


4Trprc" 


sin(u)XT) 

sin(wxL) 

sin(uxw) | 

wxT 

T 

WXW 

(6  ) 


In  (6)  Mq  is  moment,  p is  density,  r is  distance,  c is  either  P or  S 
velocity  and  is  the  radiation  pattern  (given  by  Haskell,  1964). 

XT  and  Xy  are  duration  times  associated  with  length  and  width,  respec- 
tively and  determined  by  fault  geometry  and  position  of  the  observer. 


*L  = 

|L(l/vR  - cos  0/c)/2| 

(7) 

XW  = 

| W(cos  sin  0)/(2c)J 

(8) 

Xx  = 

t/2 

(9) 

SPECTRAL  CHARACTERISTICS 

For  the  present,  let  us  adopt  (in  slightly  modified  form)  the 
similarity  relations  given  by  Kanamori  and  Anderson  (1975b). 


W 

L 


n 

L 


C^  = const 


(10) 


const 


(11) 
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er 

L 


const 


(12) 


(10)  is  the  condition  of  the  geometrical  similarity;  (11)  and  (12)  imply 
constant  stress  drop  and  constant  effective  stress. 

We  will  select  values  of  the  constants  which  seem  to  be  good 
averages  of  observational  data.  We  found  that 


L = 2W 


(13) 


seemed  to  be  the  approximate  average  of  the  empirical  data.  When  we 
substitute  L = 2W  into  (5)  and  set  6 =4.0  km/sec  we  get 

t = [16  /l2/2  / (7tt3/2  • 4)]  = . 0726L  (14) 

where  t is  in  seconds  and  L is  in  kilometers. 

We  could  use  (11J  directly  to  get  a scaling  relation  between 
fault  displacement  and  length.  In  practice  though,  most  estimates  of 
D in  Table  1 come  from  dividing  by  pS,  so  it  seems  better  to  relate 
moment  directly  to  length.  Setting  L = 2W  in  (4)  gives  moment  In  terms 
of  fault  length  and  stress  drop. 

Mq  = lAo  • 16/(7 (2tt) 3/2)  = (1.45  x 1020)  L3Ao  (15) 

where  Mq  is  in  dyne  cm,  L is  in  km  and  A o in  bars. 
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From  (6)  we  can  isolate  a spectral  factor,  dependent  only  on  fault 
parameters  and  frequency. 


sin(u>xT) 

sin(u>xL) 

sin(u>xw) 

WXT 

WXL 

<0XW 

(16) 


The  L term  follows  from  the  similarity  relation  Mq  ~ L3.  When  A(co) 
is  multiplied  by  stress  drop  and  the  constant  in  similarity  equation 
(15)  we  get  the  source  moment  rate  spectrum. 


Equation  (16)  and  the  factors  (7)-(9)  are  well  known  results  for 
the  case  of  a rectangular  fault  in  a whole  space.  These  expressions 
can  also  be  applied  directly  to  the  case  of  a rectangular  fault  with 
horizontal  rupture  propagation  in  the  earth.  The  form  of  the  expressions 
remains  identical,  but  c now  should  be  interpreted  as  the  apparent  velocity 
along  the  earth’s  surface. 

For  body  waves  c = v /sin  i,  where  v is  the  near-field  P or  S 

c c 

velocity  and  i is  the  takeoff  angle  of  the  teleseismic  ray  from  the  focal 
sphere.  This  can  be  understood  physically  by  invoking  reciprocity. 

Signals  from  a source  at  the  position  of  the  teleseismic  receiver  would 
be  picked  up  (L  cos  0)/ (2c)  sooner  at  the  end  of  the  fault  than  at 
the  center.  Thus  for  the  case  of  infinite  rupture  velocity,  this  is  the 
difference  between  arrivals  at  the  receiver  from  the  end  and  center  of 
the  fault.  This  type  of  geometrical  interpretation  can  be  applied  to 
both  (7)  and  (8),  so  that  these  factors  are  seen  to  be  the  difference 
in  arrival  times  obtained  from  geometrical  optics.  Ben-Menahem  (1962) 
gives  a more  rigorous  derivation  of  this  result. 
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For  surface  waves  (7)  is  the  well  known  directivity  factor  first 

given  by  Ben-Menahem  (1961).  If  we  neglect  the  variation  of  the  excitation 

function  with  depth,  (8)  is  the  factor  for  the  effect  of  fault  width  on 

the  surface  wave  spectrum.  In  both  (7)  and  (8),  c is  the  (frequency 

dependent)  surface  wave  phase  velocity.  The  geometrical  interpretation 

of  (7)  and  (8)  as  phase  delay  between  "arrivals"  from  the  center  of  the 

fault  and  the  ends  is  the  same  as  for  body  waves. 

Typical  values  of  c for  teleseismic  P waves  might  be  14  km/sec,  while 

for  surface  waves  4 km/sec  is  appropriate.  If  rupture  velocity, 

v = 2.5  km/sec,  then  for  surface  waves  1/c  will  be  of  the  same  order 
R 

of  magnitude  as  1/v  . As  0 varies  from  0 to  2n,  y will  range  from 
K L 

0.15L  to  0.65L.  Thus  a horizontally  propagating  rupture  will 
cause  a large  directivity  effect.  On  the  other  hand,  for  body  waves 
from  the  same  source,  cos  0 = sin  i,  so  because  of  the  relatively  steep 
takeoff  angles  of  teleseismic  rays,  it  is  reasonable  to  assume  | cos  0 | < 0 . 5 . 

This  assumption  leads  to  the  conclusion  that  y^  will  vary  only  from  0.36L 
to  0.44L.  There  will  be  only  a small  azimuthal  dependence  (i.e.  directivity 
effect)  of  the  teleseismic  body  wave  pulse.  This  implies  that  one  can 
infer  the  nature  of  a horizontal  rupture  propagation  much  more  easily 
from  surface  waves  than  teleseismic  hody  waves. 

Because  c is  the  apparent  velocity,  rather  than  the  near-field  P or 
S velocity,  it  is  inadequate  to  present  only  a single  spectrum  representing 
the  effect  of  the  seismic  source,  as  was  done,  for  example,  by  Aki  (1967,1972). 
There  is  a separate  "source  spectrum,"  A(m),  for  each  body  wave  phase  or 
surface  wave  mode,  because  of  the  different  values  of  c.  The  source 
spectrum  depends  on  azimuth,  source  dimensions,  and  the  phase  velocity,  c. 

Both  the  length  factor,  (7)  and  the  width  factor,  (8),  arc  different  for 
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each  mode.  In  the  next  section  we  will  average  and  over  all 

azimuths.  In  these  averages,  the  value  of  c will  affect  only  x,,,  even 

W 

though  both  factors  are  affected  at  nearlv  "Articular  azimuth. 

Both  (6)  and  (16)  completely  neglect  the  effect  of  the  earth's 
transfer  function  on  observed  seismic  waves.  If  one  were  to 
calculate  synthetic  seismograms  for  an  individual  earthquake,  it  would 
be  necessary  to  consider  the  earth's  response  and  the  earthquake  source 
parameters,  e.g.,  Langston  and  Helmberger  (1975)  for  body  waves  or 
Harkrider  (1964,  1970)  for  surface 

waves.  Langston  and  Helmberger  (1975)  have  demonstrated  that  sP  and 
pP  phases  play  a crucial  role  in  the  "P  wave"  from  shallow  earthquakes. 
Similarly  one  must  consider  the  surface  wave  excitation  functions  and 
the  source  mechanism  to  calculate  accurate  Rayleigh  and  Love  amplitudes. 

In  this  paper  we  consider  trends  among  events,  rather  than 
accurate  determination  of  parameters  of  particular  events.  Therefore, 
we  assume  that  the  effect  of  the  earth  structure  averages  out 
when  we  construct  scaling  relations.  Thus  we  will  use  (16)  to  get 
relations  between  m^  and  Mg,  log  L and  Mg,  log  Mq  and  M^. 
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AVERAGE  SPECTRA 


We  now  want  to  find  average  asymptotic  forms  for  log  A(w)  from 
(16).  In  particular  we  require  expressions  for  teleseismic  P phases 
(from  which  will  be  determined)  and  for  20-second  Rayleigh  waves 
(from  which  we  find  Mg).  For  both  cases  we  will  find  average  values 

of  x and  x which  take  the  direction  of  radiation  into  account.  In 

AL  u> 

making  our  approximation,  we  will  replace  |(sin  X)/X|  by  one  for 
X < 1 and  by  X ^ for  X S 1. 

Takeoff  angles  of  teleseismic  body  waves  are  nearly  vertical.  We 
will  adopt  the  approximation  that  the  rays  takeoff  straight  down. 

Thus,  for  body  waves,  we  set  6 = — in  (7)  and  (8).  Also,  in  that  case, 
| cos  4>  | = sin  6,  where  6 is  dip  angle  of  the  fault  plane.  Using  these 
values,  average  spectral  factors  for  body  waves  are 


<*L> 

body 


(17) 


and 


<XW> 


W sin  6/ (2c) 


(18) 


body 


For  surface  waves  we  will  average  and  x,,  for  9 = 0 to  0 =*  2 tt . 

L W 

On  the  earth's  surface  we  get  |cos  <j>  | = cos  6.  Thus  wo  get 


<*L> 

surf 


(19) 


and 


<XW>  = W C0S  (20) 

surf 

Comparison  of  (17)  and  (18)  with  (19)  and  (20)  shows  that  the  average 

corner  frequency  due  to  fault  length  will  be  the  same  for  body  waves 

and  surface  waves,  but  that  the  corner  frequencies  due  to  width  will  be 

different.  This  difference  affects  the  high  frequency  spectrum  only 

since  the  average  corner  frequency  for  width  is  higher  than  that  for 

rise  time  or  length.  Note  that  we  have  assumed  that  rupture  propagates 

parallel  to  the  earth's  surface  to  obtain  (17)-(20). 

Before  calculating  numerical  values  for  (17)-(20)  we  must  fix  V , 

5 and  c.  Also  we  will  use  (13)  to  relate  L to  W.  We  will  set 

VR  = 2.88  km/sec,  6 = 4.5°,  c = 14  km/sec  for  body  waves  and 

c *=  3.9  km/sec  for  surface  waves.  (For  the  earth,  c must  be  the 

appropriate  phase  velocity,  not  the  S wave  velocity.  Neglecting  the  frequency 
dependence  of  surface  wave  phase  velocity  is  a reasonable  approximation.) 

From  (9)  and  (14) 

X * 0.0363  L - C L C21) 

IT 
(17)  and  (19)  both  give 


(23) 
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For  surface  waves  we  get 

< Xw>  = 0.0289  L = C L (24) 

surf  “ 

We  now  can  approximate  the  logarithm  of  A(w)  from  (16). 
log  A(w)  = 3 log  L for  u>  < (C  L)”1 

log  A(u)  = 2 log  L - log  to  - log  CT  for  (C  L)"1  < to  < (C  L)_1 

(25) 

log  A (to)  =-  log  L - 2 log  to  - log  (ClCt)  for  (C^L)-1  < to  < (OX)"1 
log  A(to)  = -3  log  to  --  log  (C^C^)  for  (CyL)"1  < to 

The  spectra  from  these  relations  are  plotted  in  Figure  4.  Note 

_2 

that  the  body  wave  spectra  have  a much  longer  interval  of  to  decay  than 

the  surface  wave  spectra.  Calibration  of  these  curves  with  M and  M 

s o 

is  discussed  below.  Also,  note  that  to  is  in  radians  in  (25),  but  in 
hz  in  Figure  4. 

The  asymptotic  spectral  amplitudes  given  by  (25)  are  very  similar 
to  the  results  obtained  by  Kanamori  and  Anderson  (1975b).  They  used 
the  same  asymptotic  approximation  for  sin  X/ X in  conjunction  with 
Haskell’s  (1964)  spectral  expression.  Since  Haskell's  expression 
ignores  the  effect  of  width  on  the  spectrum,  the  results  of  this  paper 
differ  from  Kanamori  and  Andersons’  only  at  frequencies  above  the  corner 
frequency  for  width.  For  example,  the  model  in  this  paper  predicts 
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constant  20  sec  spectral  amplitude  for  faults  longer  than  110  kra  while 
Kanamori  and  Andersons'  predicts  amplitudes  which  increase  linearly  with 
L.  As  a result,  their  model  predicts  log  L for  large  events, 

while  this  paper  predicts  Mg  = const. 

-3 

Even  though  both  spectra  in  Figure  4 have  an  eventual  u > asymp- 

tote, they  are  quite  different  than  Aki's  (1967)  w-cube  model.  Aki's 
models,  as  a result  of  his  assumption  that  vk^  = k^,  had  only  a single 

corner  frequency.  His  w-cube  model  makes  a fairly  abrupt  transition 
0 -3 

from  w to  u behavior.  The  spectra  presented  here,  particularly  the 

0 -3 

body  wave  spectrum,  show  a gradual  transition  from  w to  w asymptotes. 
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n^-M  RELATION 

Changes  in  the  definition  of  the  body  wave  magnitude  scale  have 
resulted  in  a large  amount  of  confusion  today.  Gutenberg  and  Richter 
(1942)  extended  the  body  wave  magnitude  scale  from  local  events  to 
fairly  distant  events  which  were  recorded  on  Wood-Anderson  and  strong 
motion  torsion  instruments. 

Gutenberg  (1945)  introduced  a scheme  for  differing  only  in 
minor  details  from  the  summary  in  Richter  (1958).  He  determined  m^ 
from  the  instruments  available  in  1945,  which  were  mostly  broadband 
mechanical  types.  Gutenberg  (1945)  stated  that  "the  average  period  of 
P waves  in  teleseisms  is  about  4-6  seconds."  In  general,  Gutenberg 
did  not  publish  the  period  of  the  P waves  he  used  in  determining  m^, 
but  from  a preliminary  examination  of  his  unpublished  data  it  seems 
that  many  of  his  amplitudes  were  obtained  at  periods  of  4-10  sec. 

Gutenberg  and  Richter  (1956)  published  their  final  version  of  the 
relation  between  m^  and  Mg : 


m = 0.63  M + 2.50 
b s 

Mg  = 1.59  rr^  - 3.97 


(26a) 

(26b) 


Their  primary  reason  for  deriving  this  relation  was  to  facilitate  the 
construction  of  a "unified  magnitude  scale."  Investigators  at  that 
time  apparently  viewed  the  discrepancy  between  the  two  magnitude  scales 
as  an  experimental  error,  rather  than  a fundamental  effect 

of  the  seismic  source  spectrum.  This  view  was  not  unreasonable  at  the 
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time  because  m^  was  measured  at  periods  differing  only  by  a factor  of 

2 to  5 from  M and  modern  source  theories  had  not  yet  been  developed.  In  any 
s 

case,  Gutenberg  and  Richter  found  the  (body  wave)  magnitude  mg,  corres- 
ponding to  a given  Mg , by  using  (26a).  They  then  took  a weighted 
average  of  m^,  the  actual  body  wave  magnitude,  and  to  obtain  m,  the 
unified  magnitude.  Later  Richter  (1958)  published  values  of  unified  (surface 
wave)  magnitude,  M,  which  he  obtained  by  converting  m to  M using  (26b). 

In  retrospect,  unified  magnitude  was  inappropriate,  since  it  now  is  clear  that 
for  all  seismic  source  theories  and  represent  different  parts  of  the 
apectrum  which  are  not  related  by  a factor  independent  of  fault  length. 

m,  determinations  by  the  USCGS  (later  NOAA  and  now  the  USGS) 
differ  markedly  from  those  used  by  Gutenberg  and  Richter.  USCGS  values 
for  m^  use  (1) , but  A and  T are  measured  on  the  WWSSN  short  period 
instrument,  which  is  sharply  peaked  at  0.5  sec.  T nearly  always  is 
about  1 sec  in  WWSSN  magnitude  determinations.  Thus  WWSSN  magnitudes 
are  based  roughly  on  1 sec  spectral  amplitude.  On  the  other  hand, 

Gutenberg  and  Richter  determined  m^  for  many  events  at  about  5 sec, 
with  even  larger  T for  the  largest  events.  Therefore,  it  is  wrong  to 
take  the  Gutenberg-Richter  m^  as  being  related  to  spectral  amplitude  at 
any  one  particular  period. 

MODELING  m -M 
b s 

Akl  (1967)  proposed  two  statistical  models  of  seismic  sources, 

_2 

an  "in-square"  model  (which  decayed  as  u>  at  high  frequencies)  and  an 

-3  ' 

"u-cube"  model,  after  Haskell  (1966)  (which  decayed  as  uj  ).  Aki 
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compared  these  two  models  by  calculating  spectral  ratios  for  similar 
events  and  by  calculating  the  relation  of  to  Mg  predicted  by  each 
model. 

Aki  calculated  Mg  by  adding  a constant  to  the  logarithm  of 

spectral  amplitudes  at  20  sec.  The  constant  was  chosen  to  give  the 

best  agreement  between  theoretical  and  observational  spectral  ratios  of 

pairs  of  similar  earthquakes  studied  by  Berckhemer  (1962).  After 

fixing  the  additive  constant  for  M , Aki  then  defined  a similar  relation 

s 

for  ra^.  He  set  = const  + (.71  .83)  log  A(1  sec)  and  found  the 

constant  which  would  make  m^  * when  M = 6.75.  The  coefficient 

of  log  A(1  sec)  comes  from  a correction  for  duration. 

Aki  (1967)  calibrated  his  curves  for  the  w-square  and  co-cube 
models  in  this  way.  He  then  compared  the  m^-l^  curves  predicted  by  the 
models  to  the  Gutenberg-Richter  m^-M^  relation  (26).  He  suggested 
that  the  excellent  agreement  of  the  co-square  model  with  (26)  strongly 
supported  it,  over  the  co-cube  model.  Unfortunately,  his 
theoretical  m^-M^  curve  was  based  on  1-sec  spectral  amplitudes,  while 
(26)  was  derived  from  mostly  4 to  10  sec  data.  Actually  it  seems  that  Aki's 
support  for  the  co-square  model  was  incorrect.  The  WWSSN  m^-M  data 
(based  on  1 sec  m^) , discussed  below,  disagree  with  the  co-square 
model . 

The  approach  in  this  paper  is  to  match  rn^-M^,  log  S-M^,  log  M^-M^ 
and  spectral  ratio  data  simultaneously,  adjusting  the  two  free 
parameters  to  get  good  overall  agreement  with  the  data. 

A least  squares  solution  is  not  particularly  appropriate 
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because  of  the  large  number  of  parameters 

and  the  lack  of  similarity  (e.g.,  different  stress  drops) 

found  when  earthquakes  are  examined  in  detail. 

m^  is  approximated  by  a constant  plus  log  A(1  sec)  and  Mg  by 

another  constant  plus  log  A(20  sec).  A(w)  was  found  using  (25)  with 

the  constants  in  (21)-(24).  After  several  trials  the  additive  constants 

for  m,  and  M were  determined  to  be  C =4.30  and  C„.  =2.97.  To 

^ s “b  Ms 

obtain  seismic  moment  as  a function  of  L,  it  was  necessary  to  assign 

Ao  for  use  in  (15).  Kanamori  and  Anderson  (1975b) 

found  that  stress  drops  are  10-30  bars  for  most  interplate  earthquakes 
and  30-100  bars  for  most  intraplate  earthquakes,  so  A a = 50  bars  was  used. 

Clearly  it  is  not  exactly  correct  to  get  and  directly  from 
spectral  amplitudes.  A more  accurate  approach  would  be  com- 
puting synthetic  seismograms  and  then  measuring  m^  and  Mg  as 
it  is  done  for  data.  For  this  study,  using  spectral  amplitudes  seems 
to  be  an  acceptable  approximation.  We  plan  to  measure  m^  and 
directly  from  synthetic  seismograms  in  a future  study. 

Archambeau  (1975)  thoroughly  discusses  the  differences  between 
time  domain  and  frequency  domain  estimates  of  m^  and  Mg.  His  study 
stressed  the  very  small  differences  which  are  crucial  in  the  context  of 
seismic  discrimination.  Tn  general  though,  his  study  supports  the 
applicability  of  using  spectral  amplitudes  to  estimate  m^  of  small  and 
moderate  events.  Probably  any  discrepancy  between  spectral  amplitudes 
and  time  domain  amplitudes  is  most  severe  for  larger  events.  Tn  later 


work  we  will  calculate  the  size  of  that  effect. 
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n^-M  DATA 

Two  kinds  of  m^-M  data  are  plotted  in  Figure  5.  Points  below 
the  solid  line  midway  up  the  figure  are  from  a study  of  almost  one 
thousand  events  by  Evemden  (1975).  Each  point  is  the  average  value 
of  for  all  earthquakes  with  that  value.  Because  Evemden's 
values  are  an  average  of  0.3  lower  than  the  USGS  values,  0.3  is  added 
to  the  values  before  plotting  them.  Points  above  the  line  are 
values  for  individual  events  since  raid-1963  as  listed  in  Table  1. 

Data  shown  in  Figure  5 are  in  general  agreement  with  a study  of 
the  m^-M^  relation  by  Nagamune  (1972).  Nagamune  fitted  two  straight 
lines  to  two  years  of  WWSSN  data.  He  found 


Mg  = 1.89  m^  - 4.62 


M > 5.73 
s 


M = 1.05  m - 0.02 


M < 5.73 


The  latter  equation  comes  from  a study  of  small  events  mostly  in 
Hokkaido.  Magnitudes  in  the  latter  equation  are  very  similar  to  m^ 


and  M . 

s 


m^-M  curves  from  two  models  are  plotted  in  Figure  5.  The  curve 

on  the  right  is  the  m - M relation  predicted  by  Aki  (1972),  which  is 

b s 

based  just  on  log  A(l),  without  any  correction  for  duration.  It  can 
be  seen  that  all  the  data  lie  substantially  to  the  left  (smaller  m^) 
of  the  u-squarc  curve.  Inclusion  of  a duration  correction  for  large 
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events  would  not  affect  tlie  basic  conclusion  that  the  w-square  model 
does  not  agree  at  all  with  the  data. 

The  lefthand  curve  is  derived  from  the  Haskell  model  presented 
in  this  paper.  It  can  be  seen  that  the  predicted  m^-M^  curve  is  gen- 
erally in  good  agreement  with  the  data.  It  would  have  been  better  to 
have  averaged  the  value  of  m^  for  all  earthquakes  with  a particular 
for  all  WWSSN  events  for  sever. il  years  .rather  than  present  just  a few 
data  points. 

The  Evernden  data  have  a slope  of  one  for  events  with  im  smaller 

3 

than  bk  while  the  predicted  curve  has  a slope  of  — for  > 4.2, 

which  clearly  disagrees  with  the  data.  The  large  events  are  too 
scattered  to  warrant  a definite  conclusion,  but  the  predicted  maximum 
m^  of  6.0  is  probably  0.3  or  0.4  too  small.  This  discrepancy  may  be 
due  to  use  of  spectral  amplitudes  instead  of  time  domain  amplitudes. 
Also,  it  was  assumed  above  that  m^  was  always  based  on  1 sec  observa- 
tions, but  this  is  not  strictly  true.  The  Portuguese  earthquake  oi 
1969  (number  34  in  Table  1)  has  m^  of  7.3,  the  largest  of  any  of  the 
events  in  Table  1.  The  average  T for  this  event  was  1.77  sec;  the 
Haskell  model  predicts  that  if  had  been  determined  at  1 sec,  it 
would  be  0.5  smaller.  A systematic  variation  of  T as  a function  of 
m^  could  account  for  part  of  the  difference  between  the  theoretical 
curve  and  observations.  Another  possible  explanation  ot  the  differ- 
ence may  be  heterogeneity  of  the  source  mechanism.  This  possibility  is 


discussed  later  in  more  detail. 
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HIGH  FREQUENCY  SPECTRA 
-3 

The  Haskell  model,  which  has  w high  frequency  decay,  moment 

proportional  to  L and  "comer  frequency",  w , (for  fixed  source- 

receiver  geometry,  source  similarity  and  source  mechanism)  proportional 

to  L 1,  is  a particular  member  of  a general  class  of  models  having 

those  properties.  Following  an  argument  first  suggested  by  Savage 

(1972) , note  that  for  many  source  models  the  area  radiating  energy  to 

2 

a far-field  observer  will  appear  to  grow  as  t and  dislocation  from 

that  area  will  grow  linearly  with  t.  The  far-field  pulse,  which  is 

2 

the  time  derivative  of  the  moment  function  will  grow  as  t , giving 

_3 

(Bracewell,  1965,  p.  144)  w high  frequency  spectral  decay  (assuming 

2 

the  t onset  is  the  most  abrupt  discontinuity).  Many  models  will  also 
have  a "comer  frequency"  proportional  to  L \ where  E is  some  char- 
acteristic source  dimension  of  that  model.  Finally,  most  models  give 

3 

M L , where  L is  a source  dimension.  For  all  models  meeting  the 
o 

above  three  requirements,  high  frequency  spectral  amplitudes  will 
behave  as 

A(m)  ~ M (mr,/w)^  ~ L^(I.  Vw)^  ~ u>  ^ 
o C 

Thus,  all  events  with  fixed  geometry  and  source  type  will  share  a 

common  high  frequency  asymptote  which  is  independent  of  L.  Therefore 

the  conclusion  that  m,  , and  for  much  larger  events,  M , will  have  a 

b s 

maximum  value,  applies  to  a more,  general  class  of  models.  For  example. 
Minster  (1973),  derived  m, -M  curves  (with  a similar  shape  to  the 


curve  from  the  Haskell  model  in  Figure  5)  from  an  Archambeau  type 
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(volume)  source  (also  u falloff),  although  he  did  not  calibrate  them 
against  m^-M  data.  Minster's  results  also  predict  maximum  values  of 
m^  and  M . 

Many  investigators,  such  as  Richards  (1973),  Dahlen  (1974),  Sato 

and  Hirasawa  (1973)  and  Madariaga  (1975)  have  outlined  crack  models  for 

-3 

which  the  initial  rupture  contributes  an  u high  frequency  spectrum 

while  a "stopping  phase"  caused  by  simultaneous  cessation  of  fracture 

-2 

everywhere  on  the  fault  contributes  w and  therefore  dominates  the 

high  frequency  spectrum.  If  such  models  are  applicable,  ra^,  which  is 

based  on  the  initial  rupture,  would  still  have  a maximum  value,  but 

M would  not. 
s 

M VERSUS  FAULT  AREA 
s 

Figure  6 is  a plot  of  fault  area,  S,  (taken  from  Table  1)  and  M^. 
The  predicted  M^-log  S curve  derived  from  (21)  — (25)  agrees  quite  well 

with  the  data.  Note  that  the  theoretical  curve  has  four  different 

2 

segments.  For  small  earthquakes , up  to  = 6.76,  the  slope  is  — . 

From  Mg  = 6.76  to  = 8.12,  the  region  in  which  most  moderately  large 

earthquakes  are  clustered,  the  slope  is  1.  There  is  a small  section 

for  which  the  slope  is  2,  from  M = 8.12  to  M = 8.22.  After  M = 8.22, 

s s s 

the  largest  value  of  M for  this  calibration  of  the  Haskell  model,  the 
slope  is  infinite  (e.g.,  S increases  with  no  further  increase  in 
magnitude. ) 

There  is  a systematic  difference  between  the  interplate  (closed 
circles)  and  intraplate  (open  circles)  in  Figure  6.  Half  ot  the 


V-25 


25 

intraplate  events  fall  below  the  predicted  M - log  S curve,  while 

nearly  all  interplate  earthquakes  are  above  the  curve.  Kanamori  and 

Anderson  (1975b)  showed  that,  at  least  in  the  region  with  slope  one, 

this  meant  intraplate  events  had  a higher  apparent  stress. 

Utsu  and  Seki  (1954)  found  the  empirical  relation  log  S = 1.02  M-4.01. 

Their  I!  is  Japan  Meteorological  Agency  (JMA)  magnitude,  which  is  roughly 

2 

equivalent  to  M ,and  S is  in  km".  For  the  unit  slope  part  of  the 

M - log  S curve  in  Figure  6 (M  > 6.76)  the  Utsu-Seki  relation  predicts 
s s 

about  five  tines  the  fault  area.  This  may  be  due  to  the  way  Utsu  and 
Seki  apparently  determined  fault  area.  They  used  an  area  encompassing 
nearly  all  the  aftershocks,  rather  than  the  one-day  aftershock  zone 
which  seems  to  give  much  better  agreement  with  observed  fault  dimensions  for 
earthquakes  on  continents.  Bath  and  Duda  (1964)  proposed  the  relation 
log  S = 1.21  M - 5.05,  based  on  a study  of  six  earthquakes  from 

/ • , 2 ‘ ‘ ’ i , 

different  regions.  Bath  and  Duda  used  S as  aftershock  area  (m  km  ),  not  fault 

plane  area,  so  this  is  basically  similar  to  Utsu  and  Seki’s  result. 

Chinnery  (1969)  summarizes  a number  of  efforts  to  find  a single  linear 

relation  between  M and  the  logarithm  of  other  fault  parameters, 
s 

M VERSUS  MOMENT 
s 

The  data  of  Kanamori  and  Anderson  (1975b)  show  that  Act  = 50  bars 

is  a good  average,  about  halfway  between  values  for  interplate  and 

intraplate  events.  Using  Aa  = 50  and  (15),  we  find  moment (in  dyne  cm) 

, 21  3 

is  related  to  fault  length  (in  km)  by  M ~ (7.26  * 10  )L  or 

log  M = 21.9  + 3 log  L.  It  was  shown  above  that  M n log  L,  where,  n 
o s 
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varies  between  0 and  3,  as  can  be  seen  from  the  surface  wave  spectra 

in  Figure  4.  Therefore,  for  small  earthquakes  the  M^ilog  slope  is 

one;  for  very  large  events  (Mg  constant)  the  slope  is  infinite 

(M  increases  but  M is  already  at  a maximum.), 
o s 

The  log  M - data  from  Table  1 are  plotted  in  Figure  7.  Most 

of  the  moderate  sized  events  (Mg  from  6.76  to  8.11)  fall  on  the  slope 

3 

1.5  portion  of  the  curve  (log  ~ — M^).  This  part  of  the  curve 

corresponds  to  cases  where  20  sec  spectral  amplitudes  are  measured  on 

the  u>  ^ part  of  the  spectrum.  Because  the  corner  frequency  for  width 

is  only  slightly  greater  than  that  for  rise  time,  the  slope  3 

(Mq  ~ 3 log  Mg)  region  is  very  small,  extending  only  from  Mg  = 8.12  to 

M = 8.22.  Beyond  that,  slope  is  infinite.  The  data  agree  quite  well 

with  the  theoretical  curve.  As  in  Figure  6,  intraplate  events  tend  to 

have  smaller  Mq  for  a given  M , corresponding  to  higher  apparent  stress. 

Aki  (1972)  showed  that  his  w-square  model  also  agreed  well  with 

vs.  log  data.  Brune  and  co-workers  (Brune  and  King,  1967; 

Brune,  1968;  and  Brune  and  Engen,  1969)  presented  a magnitude  scale 

based  on  100  sec  surface  wave  amplitude.  They  then  assumed  log 

log  A(100)  and  fit  two  segments,  each  with  slope  1,  to  the  data. 

Because  of  their  different  definition  of  M , their  results  cannot  be 

o 

directly  compared  to  this  paper. 

Data  in  this  section  show  that  when  M is  larger  than  about 

o 

2 ^ 

10  dyne  cm,  M reaches  its  maximum  value.  It  is  important  to  con- 
s 

sider  this  when  discussing  the  "maximum  credible  earthquake"  likely  to 
occur  in  a particular  area.  The  earthquake  size  may  be  specified  in 
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terms  of  M for  most  earthquakes,  but  when  the  moment  approaches  10  , 

magnitude  no  longer  is  a valid  parameter  for  specifying  earthquake 
size.  Whenever  the  maximum  crediole  earthquake  is  in  this  range,  e.g., 
as  is  probably  the  case  in  discussing  the  Alaskan  pipeline,  moment, 
not  Mg,  should  be  the  parameter  used. 


SPECTRAL  RATIOS  OF  SIMILAR  EVENTS 


Berckheraer  (1962)  studied  spectral  ratios  of  earthquakes  with 
roughly  the  same  location  and  mechanism.  In  theory,  the  spectral  ratio 
method  eliminates  the  effect  of  earth  structure  and  leaves  only  effects 
due  to  the  difference  in  source  spectra.  Aki  (1967)  used  BerckhemerTs 

data  to  determine  the  relation  between  M and  comer  period  for  the 

8 

w-square  and  m-cube  models. 

Berckheraer * s original  data  and  Aki’s  theoretical  curves  are  shown 
in  Figure  8 together  with  the  theoretical  curve  from  the  model  in  this 
paper.  Both  models  seem  to  agree  fairly  well  with  the  data.  Perhaps 
Aki's  fits  slightly  better.  Berckhemer  presented  six  pairs  of  spectra, 
of  which  only  four  are  presented  here.  The  remaining  two  pairs  used 
smaller  earthquakes,  involving  mostly  short-period  data,  which  probably 
are  less  reliable.  No  attempt  at  fitting  these  two  pairs  was  made. 

Tsujiura  (1973)  published  spectral  ratio  data  for  many  pairs  of 
earthquakes.  Most  of  his  events  could  be  fit  by  both  Aki's  w-square 
model  or  Aki's  (1972)  version  of  Brune's  "w-model",  although  usually 
one  model  or  the  other  fit  somewtiat  better.  There  were,  however,  two 
pairs  of  events  from  the.  Aleutians  which  had  spectral  ratios  that  were 
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unusually  flat  and  could  not  be  fit  by  either  model.  Tsujiura's 
spectral  ratio  data  have  not  yet  been  compared  to  the  model  in  this 
paper . 


DISCUSSION 

The  Haskell  model  with  parameters  (21)— (24 ) is  in  general  agree- 
ment with  m^-M^  data  (Fig.  5),  - log  S data  (Fig.  6),  Mg  - log 

data  (Fig.  7)  and  spectral  ratio  data  (Fig.  8).  The  most  serious 
discrepancy  between  the  data  and  the  model  comes  in  Figure  5.  On  one 
hand,  the  maximum  value  of  m^  is  probably  several  tenths  too  small. 

On  the  other,  the  data  seem  to  have  a slope  of  about  one  up  to  = 5-i, 
while  the  curve  from  the  model  has  slope  one  only  up  to  = 4.19. 

This  phenomenon  could  be  explained  if  most  earthquakes  are  complex 
sources  with  the  first  burst  of  energy  coining  from  a smaller,  sub- 
stantially higher  stress  drop  source,  than  the  average  of  the  whole 
earthquake.  If  this  is  the  case,  then  m^  would  be  measured  on  a flat 
or  flatter  part  of  the  spectrum  than  one  would  expect  for  the  earth- 
quake as  a whole. 

Burdick  and  Mellman  (1976)  have  suggested  that  for  the  Borrego 
earthquake  of  1968  most  of  the  body  wave  energy  came  from  a source 

region  with  radius  of  8 km,  giving  about  half  the  area  shown  in  Table  1. 

2 7 

Since  they  also  found  a higher  moment,  0.112  * 10"  dyne  cm,  their 
stress  drop,  96  bars,  is  about  4 times  the  value  in  Table  1,  taken  from 
Hanks  and  Wyss  (1972).  Tucker  and  Bruno  (1975)  also  suggested  that 
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sources  showed  a smaller  high  stress  drop  event  superimposed  on  the 
overall  average  event.  The  m^-M^  data  in  Figure  5 agree  with  the 
possibility  of  the  initial  fracture  having  higher  stress  drop  than  the 
bulk  event,  but  certainly  do  not  prove  that  this  happens.  Other 
explanations,  e.g. , the  whole  model  is  wrong,  or  similarity  is  wrong 
for  smaller  shocks,  are  equally  admissible. 

SUMMARY 

The  following  scaling  relations  relating  width  and  rise  time  to 
length  and  fault  area  have  been  given; 

h = 2W 

r - 16S1/2  / (7tt3/23)  . 

The  relation  for  rise  time  was  derived  from  the  assumption  that  static 
stress  drop  and  dynamic  effective  stress  are  equal;  agreement  of 
theoretical  rise  times  with  the  data  tends  to  support  that  assumption. 

Averages  of  observed  rupture  velocities  show  that  = 0.723. 

The  Haskell  model  predicts  that  magnitude  will  reach  an  upper 
limit  regardless  of  further  increases  in  fault  length  and  seismic 
moment.  Moment,  rather  than  magnitude  should  be  used  to  discuss  the 
possible  size  of  great  earthquakes. 

The  "source  spectrum”  from  any  source  model  is  a function  of 
apparent  (phase)  velocity  of  the  mode  or  phase  being  considered,  as  well 
as  of  azimuth  and  source  parameters.  It  is  incorrect  to  speak  of  a 
single  "source  spectrum." 
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Theoretical  relations  between  m_  £ind  M from  the  Haskell  model 

1)  s 


are : 


= Mg  + 1.33 


M <2.86 
s 


m^  = 3^  + 2.28 


2.86  < M < A. 90 
s 


mb  = + 3-91 


4.90  < M < 6.27 

s 


m,  = 6.00 
□ 


6.27  < M 

s 


M and  fault  area  (in  km  ) are  related  by 
s 


log  S = | Mg  - 2.28 


M <6.76 
s 


log  S = M - 4.53 
s 


6.76  < M < 8.12 
s 


log  S = 2M  - 12.65 
s 


8.12  < M < 8.22 
s 


M = 8.22 

s 


S > 6080  km 


if  L = 2W  is  used. 


1 
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If  we  assume  a stress  drop  of  50  bars,  then  log  (in  dyne  cm) 

and  M are  related  by 
s 


log  M = M + 18.89 
o s 

log  M = -qM  + 15.51 
3 o 2 s 

log  M = 3M  + 3.33 
& o s 


M <6.76 
s 

6.76  < M < 8.12 

s 

8.12  < M < 8.22 
s 


M = 8.22  log  M > 28 

B O 

These  scaling  relations  fit  observed  data  quite  well.  They  should  not 
be  used  to  determine  the  value  of  a parameter  for  any  individual  earth- 
quake, since  these  "averages",  and  the  assumptions  made  to  derive  them, 
are  not  exactly  correct  for  any  single  event. 

A review  of  work  by  Gutenberg  and  Richter  reveals  that  their 
m^-M^  relation  was  derived  from  data  at  mostly  5 or  10  second  period. 
Models  such  as  Aki's  (1967)  w-square  model  which  fit  theoretical  1 sec 
data  to  the  Gutenberg-Richter  relation  are  probably  in  error. 
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FIGURE  CAPTIONS 

Plot  of  fault  length  (along  strike)  versus  fault  width  (along 
dip)  for  earthquakes  in  Table  1.  Open  circles  are  intra- 
plate events;  closed  circles  are  interplate  events.  Numbers 
refer  to  Table  1.  These  conventions  are  used  for  all  plots 
of  earthquake  data. 

Plot  of  observed  rise  time  versus  theoretical  rise  times 
from  (5). 

The  Haskell  (whole  space)  fault  model. 

Source  spectra  of  surface  waves  (on  left)  and  body  waves. 

_2 

Both  are  identical  at  frequencies  below  the  w corner 

frequency.  The  body  waves  have  a higher  width  corner 

frequency  than  surface  waves,  which  follows  from  (18),  (20), 

(23)  and  (24).  This  difference  occurs  because  teleseismic 

P waves,  which  takeoff  essentially  straight  down,  have  a 

much  higher  apparent  velocity  (phase  velocity)  than  surface 

waves.  Therefore  the  separation  between  rise  time  and  width 

-2 

corner  frequencies  (the  w part  of  the  spectrum)  is  much 
greater  for  body  waves  than  surface  waves. 

USCGS  m^  versus  M^.  Lower  data  points  are  averages  from 
Evernden  (1975),  corrected  by  adding  0.3  to  m^.  Upper  points 
(above  horizontal  line)  are  individual  earthquakes  from 
Table  1.  Dashed  curve  is  m^-M^  relation  from  w-square 
model.  Solid  curve  is  from  (21)  to  (25)  as  described  in 
the  text. 

M versus  log  S data  with  theoretical  curve  from  the 


model ’n  this  paper. 
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Figure  7. 


Figure  8. 


M versus  log  M data  with  theoretical  curve  from  t^e 
s o 

model  in  this  paper. 

Spectral  ratios  of  similar  earthquakes  from  Berckhemer 
(1962).  Numbers  above  each  figure  are  magnitudes  of  the 
larger  and  smaller  of  the  events.  Dashed  line  is  w-square 
model  and  solid  line  is  model  from  this  paper. 
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